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SECTION  I 


INTRODUCTION 


The  SCEPTRE  Circuit  Analysis  Program  has  been  developed  and 
Improved  over  a  period  of  years  on  a  number  of  contracts  led  by  the 
Air  Force  Weapons  Laboratory.  The  users'  information  for  the  1973 
version  of  this  program  is  contained  in  two  volumes,  of  which  this  is 
the  second. 

Volume  I  contains  basic  instructions  for  use  of  the  program. 
Volume  II  is  intended  to  give  a  fuller  understanding  of  the  SCEPTRE 
internal  formulation  so  that  the  more  sophisticated  user  may  take  full 
advantage  of  the  program's  capabilities  and  flexibility.  This  volume 
supercedes  AFWL-TR-72-77  for  use  of  SCEPTRE  on  S/360  equipment,  and 
retains  the  Information  from  AFWL-TR-72-77  pertinent  to  7090/94  equip¬ 
ment. 

The  1973  version  of  SCEPTRE,  Implemented  for  S/360  use,  contains 
six  additions  not  available  on  previous  versions.  There  are  four  new 
DC  options,  an  AC  analysis  capability,  and  a  convolution  capability. 

This  volume  covers  the  mathematical  formulations  for  all  of  the  previous 
SCEPTRE  features,  plus  five  of  the  six  new  ones.  The  convolution  feature 
is  covered  entirely  in  Volume  I. 


SECTION  II 


FORMULATION  AND  THEORY 


2.1  INTRODUCTION 


Any  automatic  transient  analysis  program  is  designed  to  relieve  the  user 
of  the  necessity  of  writing  and  programming  the  differential  and  algebraic 
equations  that  describe  networks.  PREDICT  and  other  programs  already  perform 
this  basic  task,  but  the  degree  of  flexibility  permitted  the  user  varies  widely 
among  programs.  SCEPTRE,  written  as  a  successor  to  PREDICT,  incorporates 
many  Improvements. 

This  section  presents  the  formulation  and  theory  that  serve  as  the  basis 
for  SCEPTRE.  Therefore,  the  discussion  is  mathematically  oriented.  Those 
features  of  SCEPTRE  that  are  not  mathematical  are  not  included  here.  For 
example,  the  extremely  useful  features  of  Rerun  and  Model  Storage  are  not  de¬ 
scribed. 

2 . 2  GENERAL  SOLUTION 


SCEPTRE  consists  of  three  separate  formulations  that  combine  to  produce 
the  general  solutions  of  a  given  network.  One  is  referred  to  as  the  DC  pro¬ 
gram.  It  has  five  options,  discussed  in  subsection  2.4,  Each  will  determine 
the  network  voltages  that  prevail  before  any  time-varying  forcing  function  is 
applied.  This  program  does  not  treat  time  as  an  independent  variable;  instead 
it  holds  time  constant,  and  iterates  on  selected  voltages.  The  output  of  a 
DC  program  may  be  obtained  independently,  or  it  may  be  automatically  used  as 
the  starting  point  for  the  transient  or  AC  program.  Thus,  the  output  of  the 
DC  program  effectively  supplies  the  initial  conditions  for  the  system  of  dif¬ 
ferential  equations  that  are  solved  by  the  transient  or  AC  program. 

The  second  program  is  called  the  transient  program.  This  program  uses 
time  as  an  independent  variable  and  solves  systems  of  differential  equations 
as  functions  of  time.  The  output  of  this  program  represents  the  transient 
response  of  a  given  network.  As  implied  above,  the  transient  program  may  be 
used  in  conjunction  with  the  DC  program,  or  it  may  be  used  by  itself  if  the 
initial  conditions  of  the  network  are  known. 

The  third  program  is  called  the  AC  program.  This  program  uses  frequency 
as  an  independent  variable  and  solves  algebraic  equations  as  a  function  of 
frequency.  The  output  of  this  program  represents  the  frequency  response  of  a 
given  network.  The  AC  program  may  also  be  used  in  conjunction  with  the  DC 
program  or  by  itself.  It  is  discussed  in  Section  2.5 

The  general  solution  procedure  described  here  concerns  the  definition  of 
the  terms,  matrices  and  procedures  that  are  common  to  all  programs.  Other 
parts  of  this  volume  will  provide  the  detailed  explanations  and  derivations. 


-  2  - 


The  first  step  in  either  program  is  the  construction  of  a  tree  (Ref.  1)  ac¬ 
cording  to  prescribed  rules,  which  differ  for  the  three  programs.  This  permits 
formation  of  a  B  matrix  that  effectively  expresses  link  voltages  in  terms 
of  tree  branch  voltages  and  tree  branch  currents  in  terms  of  link  currents. 
Figure  1  shows  a  composite  B  matrix  that  contains  all  possible  element  classi¬ 
fications  and  submatrices.  This  matrix  is  derived  in  appendix  B  •  The  element 
classifications  are  given  in  table  I. 
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Figure  1.  Composite  B  Matrix  in  SCEPTRE 


*  A  tree  is  defined  as  any  connected  network  subgraph  that  contains  all  nodes 
of  the  network  but  no  complete  loops.  All  circuit  elements  that  are  members 
of  the  tree  are  termed  tree  branches.  All  circuit  elements  excluded  from  the 
tree  are  termed  links.  A  "C"  tree  is  defined  in  this  report  as  one  in  which 
tree  members  are  chosen  in  the  preference  order  E,  C,  R  and  L.  All  current 
sources  (J)  must  be  excluded  from  the  "C"  tree.  Therefore,  these  sources 
are  links.  A  cut  set  is  defined  as  that  group  of  elements  that  would  isolate 
two  groups  of  nodes  when  removed  from  a  network. 


TABLE  I 


ELEMENT  CLASSIFICATIONS  IN  SCEPTRE 


ELEMENT 


Capacitor  links 
Resistor  links 
Inductor  links 
Capacitor  tree  branches 
Resistor  :ree  branches 
Inductor  tree  branches 
Voltage  sources 
Independent  current  sources 

Primary  current  sources  (dependent  on  voltage  across  terminals) . 

This  class  will  appear  only  in  the  derivation  of  the  initial  conditions 
and  AC  programs. 

Secondary  current  sources  (dependent  on  other  current  sources) .  This 
class  will  appear  only  in  the  derivation  of  the  initial  conditions 
program. 

Voltage  sources  that  are  dependent  on  resistor  voltages.  This  class 
will  appear  only  in  the  derivation  ol  the  transient  and  AC  programs. 

Current  sources  that  are  dependent  on  resistor  currents.  This  class 
will  appear  only  in  the  derivation  of  the  transient  and  AC  programs. 

Voltage  sources  that  are  generated  in  adjoint  calculations  due  to 
the  presence  of  secondary  current  sources  in  the  network. 


The  following  matrices  and  vectors  may  also  be  defined  as  a  result  of 
the  element  classification  in  table  I: 

1&22  -  a  diagonal  matrix  composed  only  of  resistor  links 


-  a  diagonal  matrix  composed  r.nly  of  resistor  tree  branches 


-  a  diagonal  matrix  composed  only  of  capacitor  links 


-  a  diagonal  matrix  composed  only  of  capacitor  tree  branches 

s  -  c  -1 

S44  44 

L.J2  “  A  matrix  composed  only  of  inductor  links  and  the  mutual 
inductance  between  Inductor  links 

-  a  matrix  composed  only  of  the  mutual  inductance  between 
inductor  links  and  inductor  branches 

L,fi  -  a  matrix  composed  only  of  inductor  tree  branches  and  the 
mutual  inductance  between  inductor  tree  branches 

Ggg  -  a  diagonal  matrix  composed  only  of  the  voltage  derivatives 
of  primary  current  sources 

1^,  -  vectors  composed  only  of  the  currents  or  voltages  as&u- 

ciated  with  capacitor  links 

^2*  ^2  ~  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  resistor  links 

I.j,  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 
ciated  with  inductor  links 

1^,  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  capacitor  tree  branches 

1^,  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  resistor  tree  branches 


Ig,  Vg  -  Vectors  composed  only  of  the  currents  or  voltages  asso¬ 
ciated  with  Inductor  tree  branches 

Jg,  Vg  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 
ciated  with  Independent  current  sources 

J^,  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 

ciated  with  primary  current  sources 

Jq,  Vq  -  vectors  composed  only  of  the  currents  or  voltages  asso¬ 
ciated  with  secondary  current  sources 

-  A  vector  composed  only  of  voltage  sources 

2.3  TRANSIENT  SOLUTION 
2.3.1  MATRIX  OPERATIONS 

The  state  variables  of  any  system  can  be  defined  as  the  minimum  set  of 
quantities  that  will  suffice  to  determine  all  other  quantities  in  the  system 
at  any  instant.  It  can  be  shown  that  the  knowledge  of  the  set  of  all  capaci¬ 
tor  tree  branch  voltage  (V^)  and  inductor  tree  link  currents  (I_)  is  sufficient 
to  determine  all  other  element  currents  and  voltages  and,  furthermore,  that 
this  selection  of  state  variables  will  allow  network  formulation  in  terms  of 
first-order  differential  equations.  The  starting  point  of  the  derivation  than 
is  that  quantities  V.  and  I.  are  known  and  the  list  of  unknowns  is  made  up  of 
V  ,  V2,  V  ,  V  V  ,  I  ,  1,1,  I  ,  I  ,  I  and  Vg.  In  addition,  the  deriva¬ 
tives  of  che  state  variables,  V^  And  l_,  must  be  obtained  in  preparation  for 
the  numerical  integration  routine  whicn  produces  the  updated  state  variables 
that  are  valid  at  the  next  time  increment . 

Some  of  the  equations  needed  to  solve  the  unknown  quantities  may  be  ob¬ 
tained  from  the  transient  solution  B  matrix  (figure  2).  The  B  matrix  itself 
arises  from  a  "C"  tree,  which  is  formed  by  an  E,  C,  R,  L  preference  order. 

Note  that  the  B  matrix  differs  fjom  the  composite  matrix  of  figure  1  in  that 
some  submatrices  are  zero  valued  and  that  no  distinction  is  made  between 
types  of  sources. 


For  example,  B.  _  must  be  zero  since  the  "C"  tree  preference  prohibits  the 
possibility  of  a  capacitor  link  closing  a  loop  that  contains  a  resistor  tree 
branch . 
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Figure  2.  Transient-Solution  B  Matrix  in  SCEPTRE 


form: 


-»14V4  -  B17E7 


v„  - 


-B24V4  '  625V5  "  B27E7 


-B34V4  -  B35V5  '  B36V6  '  B37E7 


'8 


'B84V4  “  B85V5  "  B86V6  “  B87E7 


(1 

(2 

(3 


V  iwupnn  in  terms  of  link  currents,  there 
Since  tree  branch  currents  can  be  written  in 

arises : 


I*  “  B14T  *1  +  B24T  X2  +  B34T  X3  +  *84*  J8 


*5  " 


b25T  l2  +  B35T  *3  +  B85  J8 


X6  - 


B36T  h  +  B86T  J8 


B17T  h  +  B27T  X2  +  B37T  X3  +  B87T  J8 


where  the  superscript  T  is  us.d  to  indicate  the  transpose  of  a  matrix. 


Two  additional  equations  may  be  obtained  £rom  differentiating  Equations  (1) 
and  (7)  yielding: 


V1  -  -BuVB17E7  (9) 

•  m  *  T  * 

*6  ■  B36  ,3+,M  J8  (l0) 

•  • 

Note  that  source  derivations  E7  and  J„  have  been  introduced  in  the  last  two 
equations*.  These  equations,  together  with  a  few  fundamental  relations,  will 
be  used  to  derive  all  of  the  network  currents  and  voltages  in  terms  of  known 
quantities. 

2. 3. 1.1  Solution  of  Resistive  Quantities 


The  resistive  quantities  of  interest  are  I^,  V^,  and  V^.  The 
fundamental  voltage  current  relations  for  resistors  permit: 

V2  '  R22I2  (ll> 

V5  '  R5  5 1 5  (12> 


may  be  explicitly  solved  for  by  manipulation  of  equations  (2),  (6),  (11) 
and  (12)  to  get 

*2  '  V1  !  'B24\  *  B27E7  -  B25R55  I  “as"  J3  +  B85T  J8|  !  <13) 


where 


«R  '  R22  +  B25R55B25  <U) 

The  significance  of  equation  (13)  is  that  the  vector  of  resistor  link  currents 
has  been  solved  in  terms  of  all  known  quantities,  since  the  right  side  of  the 
equation  is  composed  entirely  of  state  variables  and  I^»  known  sources  E^ 
and  J8  and  known  incidence  submatrices.  Once  12  is  known,  the  vectors  I5, 

V2,  and  V,.  can  be  determined  from  Equations  (6),  (11),  and  (12),  respectively. 


See  subsection  2.3.3  for  a  discussion  on  source  derivatives. 
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An  alternate  approach  may  sometimes  be  preferable.  Equations  (2),  (6), 
(11),  and  (12)  may  also  be  manipulated  to  obtain 


M, 


-1 


-Boe  G„„  B„,V.  +  B„,EJ  +  B  *  I,  +  B 


25  22  24  4  '  u2Tl 


'35  3  85 


(15) 


where 


MG  "  G55  +  B25  G22B25 


(16) 


Equation  (15)  gives  the  vector  of  resistor  tree  branch  voltages  in  terms  of 
quantities  that  are  all  known.  Then,  V.,  I.  and  I,  can  be  solved  by  Equa¬ 
tion  (2),  yielding 


and 


G22V2 


*5  “  G55V5 


(17) 

(18) 


respectively.  The  two  approaches  differ  in  the  size  of  the  matrix  to  be  in¬ 
verted.  The  first  approach  requires  the  inversion  of  a  matrix  (M^)  containing 
the  number  of  rows  and  columns  equal  to  the  number  of  class-2  elements  in  ..he 
network.  The  second  approach  requires  the  inversion  of  a  matrix  (M^)  con¬ 
taining  the  number  of  rows  and  columns  equal  to  the  number  of  class-5  elements 
in  the  network.  Networks  containing  resistors  that  are  all  constant  require 
only  one  matrix  inversion.  There  is  no  practical  difference  between  the  two 
approaches.  If,  however,  the  network  contains  at  least  one  variable  resistor, 
a  matrix  must  be  inverted  at  each  solution  time  increment.  Hundreds,  or  even 
thousands,  of  matrix  inversions  are  necessary  and  the  size  of  the  matrix  be¬ 
comes  of  significant  importance.  SCEPTRE  will  automatically  determine  which 
of  the  two  approaches  should  be  taken  for  each  individual  network. 


2.3. 1.2  Solution  of  Capacitor  Quantities 

The  capacitor  quantities  that  must  be  solved  at  each  time  step  are 
*1’  *4’  ^1  an<*  ^4'  ^ector  V,  itself  will  have  been  updated  by  the  integration 
routine  and,  hence,  will  be  known.  The  fundamental  relationships  for  capaci¬ 
tors  permit 


V4  - 

S44X4 

(19) 

V1  “ 

S11I1 

(20) 

Equations  (5), 

(9),  (19),  and 
-1  1 

(2)  may  be  combined  to  obtain 

T  T  T  I  „  '  ) 

X1  “ 

where 

MS  J-B14S44 

B24  I2  +  B34  I3  +  B84  Js|  '  B17E7  [ 

(21) 

M  - 

s 

S11  *  B14S44BUT 

(22) 
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Ac  this  poinC,  vector  I.  has  been  isolated  In  terms  of  all  known  quantities. 
There  remains  to  obtain  1^,  V^,  and  from  Equations  (5),  (1),  and  (19), 
respectively. 


2. 3. 1.3  Solution  of  Inductive  Quantities 

The  inductive  quantities  that  must  be  solved  at  each  time  step  are 

• 

I_  v,,  V,,  and  I,.  Vector  I-  will  have  been  updated  at  the  start  of  each 
tike  step  and  will  be  known,  xhe  fundamental  relations  for  inductors  permit 


V 

V 


3 

6 


'33*3  +  L36I6 

(23) 

•  • 

’63T3  +  L66T6 

(24) 

Equations  (3),  (10),  (21) ,  and  (24)  may  be  combined  to  obtain 


I-.  - 


-1 


»L  *i  -V4  -  B35V5  '  B37E7  ‘ 


B36L66B86T  4  h6B86T 


8  (25) 


where 


"l  ■  L33  +  L36B36T  +  B36L63  +  B36L66B36T 


(26) 


Now  is  written  in  terms  of  known  quantities.  Following  this,  V&  and 
1^  may  be  ontained  from  Equations  (10)  and  (23),  (10)and  (24),  and  (7),  re¬ 
spectively. 


2.3. 1.4  Solution  of  Voltage  Source  Currents  and  Current  Source  Voltages 

A  complete  list  of  possible  outputs  of  a  network  would  include  the 
current  through  voltage  sources  and  the  voltage  across  current  sources. 
These  can  be  obtained  directly  from  Equations  (8)  and  (4),  respectively, 
since  the  right  sides  of  these  equations  are  known  at  this  stage  of  the 
computational  sequence.  These  steps  complete  the  formal  derivation  of  all 
network  currents  and  voltages. 


2.3.2  SCANNING  PROCEDURE 

2. 3. 2.1  Ideal  Operation  (Submatrices  B^,  and  ■  0) 

The  series  of  matrix  operations  described  in  subsection  2.3.1  could 
very  well  be  programmed  as  they  are  to  produce  the  solution  of  the  general 
transient  analysis  problem.  All  of  the  matrix  multiplication,  addition,  etc. 
could  be  performed  at  each  time  step  to  generate  the  necessary  currents  and 
voltages.  However,  a  very  significant  improvement  in  computer  running  time 
could  be  achieved  by  a  more  efficient  utilization  of  the  information  con¬ 
tained  in  the  B  matrix. 
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Various  network  voltages  and  currents  could  be  re|d  or  "scanned"  directly 
from  the  rows  and  columns  of  the  B  matrix  respectively  .  This  can  be  put 
more  precisely  by 

\  •  -»*tb  <27) 

and 

IIB-  BTIl  (28) 


where  I^b*  II*  ^TB  and  Vl  represent  the  vectors  of  tree  branch  currents,  link 
currents,  tree  branch  voltages,  and  link  voltages,  respectively.  If  the 
vectors  and  the  B  matrix  are  partitioned  according  to  the  form  of  figure  2, 
Equations  (27)  and  (28)  lead  to  the  first  six  equations  of  subsection  2.3.1 
Quantities  V~  and  I.  are  explicitly  written  in  terms  of  known  quantities  of 
Equations  (2;  and  (6)  if  submatrix  B25  ■  0.  Once  the  resistive  quantities 
are  determined,  equation  (5)  presents  1^  in  terms  of  known  quantities  if  sub¬ 
matrix  B^  "0.  In  the  same  manner,  equation  (3)  presents  in  terms  of 
known  quantities  if  submatrix  B_,  -  0.  Clearly,  then,  the  condition  "B-,  , 
B25’  B36  "  0"  permits  solution  of  V_ ,  1^,  1^,  V  ,  I,,  and  V^  directly  by 
scanning  and  without  any  matrix  manipulation.  Quantities  I_,  V^,  I.  ,  V,, 
and  I_  c«?n  then  be  determined  from  the  operations  given  or  Implied  In  tne 
non  zero  portions  of  Equations  (11),  (12),  (21),  (24),  (19),  and  (25),  re¬ 
spectively.  Once  all  of  these  quantities  are  explicitly  obtained  in  terms  of 
known  quantities,  the  equations  are  stored,  compiled,  and  executed  at  each 
time  increment  without  recourse  to  repeated  matrix  manipulation. 


2. 3. 2. 2  Operation  When  B2J  J*  0 


In  most  large  networks,  submatrix  Bj^  does  not  equal  zero.  When  this 
happens,  quantities  V.  and  I,  cannot  be  "scanned"  out  and  either  Equation  (13) 
or  (15)  must  be  solved.  Both  of  these  equations  require  the  inversion  of  a 
matrix;  thus,  some  matrix  manipulation  must  be  done.  To  illustrate  the  pro¬ 
cedure,  assume  that  some  hypothetical  network  has  given  rise  to  the  B  matrix 
shown  in  figure  3.  The  network  is  fairly  typical  in  that  B^,  B^  ■  0,  but 
B2J  T*  0.  The  nature  of  submatrix  82^  (outlined  in  figure  3)  is  such  that 
resistors  and  contribute  no  rows  or  columns,  and  the  scanning  process 
immediately  yields 

V  _  ■  E.  -  V  ,  and  I  -  I  .  from  which  follow 

RJ  1  \jL  KD  LX  y 

rR3  ■  «1  -  VC2>/R3  and  ■  Vu 

These  quantities  will  be  updated  at  each  time  step  without  matrix  manipulation. 
The  rest  cf  the  resistive  quantities  in  the  network  may  be  solved  by  the 
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matrix  manipulation  implied  in  Equations  (13)  and  (14).  Once  the  resistive 
quantities  are  determined,  the  remaining  network  quantities  may  be  scanned 
out  and  stored.  The  only  repeated  matrix  manipulation  required  will  be  for 
the  three  resistor  link  currents  (I_  ,  I  ,  Ip  )  that  cannot  be  scanned  be¬ 
cause  B25  +  0.  K1  2  4 
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Figure  3.  B  Matrix  from  Hypothetical  Network 


2.3.3  SEMI-AUTOMATIC  SOURCE  DERIVATIVES 

The  need  for  source  time  4erivatives  in  transient  analyses  is  established 
in  Equations  (9)  and  (10)  where  and  T.  are  required.  In  situations  where 
non-zero  source  time  derivatives  ax.’  needed,  the  user  must  supply  them.  These 
situations  are  subsequently  described. 

2. 3. 3.1  Voltage  Source  Is  Variable  and  +  0 


If  B17  +  0  and  the  voltage  source  is  constant,  SCEPTRE  will  auto¬ 

matically  supply  a  zero  derivative.  In  the  case  where  a  non-zero  derivative 
is  required  and  the  user  fails  to  supply  it,  the  run  will  be  terminated  with 
a  diagnostic  message.  The  situation  is  best  recognized  by  the  presence  of 
any  circuit  loop  composed  solely  of  capacitors  and  at  least  one  variable 
voltage  source. 
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2. 3. 3. 2  Current  Source  Is  Variable  and  Bg^  i  0 


If  Bg.  t4  0  and  the  current  source  is  constant,  SCEPTRE  will  auto¬ 
matically  supply  a  zero  derivative.  In  the  case  where  a  non-zero  derivative 
is  required  and  the  user  fails  to  supply  it,  the  run  will  be  terminated  with 
a  diagnostic  message.  The  situation  is  best  recognized  by  the  presence  of 
any  circuit  cut  set  composed  solely  of  inductors  and  at  least  one  variable 
current  source. 

2.3.4  LINEARLY  DEPENDENT  SOURCES 

The  ability  of  SCEPTRE  to  correctly  process  linearly  dependent  sources 
permits  the  user  to  define  current  and  voltage  sources  that  are  linear  func¬ 
tions  of  resistor  currents  and  voltages.  This  feature  is  expected  to  be  most 
useful  for,  but  not  limited  to,  applications  involving  families  of  small- 
signal  transistor  equivalent  circuits  such  as  shown  in  figure  4. 

These  linearly  dependent  sources  require  special  treatment  because  their 
magnitudes  are  direct  functions  of  quantities  that  are  not  state  variables. 
Unless  these  sources  are  specifically  processed,  they  will  be  updated  at  the 
n^  time  step  according  to  the  values  of  their  independent  variables  at  the 
(n  -  1)  time  step  (as  they  would  be  in  PREDICT,  for  example).  This  results 
in  a  "computational  delay",  which  can  lead  to  large  errors  throughout  the 
entire  network. 

Linearly  dependent  sources  are  provided  for  in  the  general  B  matrix 
by  the  Y  and  X  classification.  When  these  sources  exist,  SCEPTRE  will  set 
up  the  B  matrix  as  shown  in  figure  5.  Under  these  circumstances.  Equations 
(2)  and  (6)  can  be  extended  to: 

V2  ■  -B24V4  '  B25V5  ‘  B27E7  '  V*  (29) 

T5  ■  B25T  J2  +  B35T  X3  +  B85T  J8  +  Bx5T  JX  (30> 

Since  any  resistor-voltage-dependent  voltage  source  must  depend  on  a  resistor 
tree  branch  voltage  or  a  resistor  link  voltage,  there  arises 

EY  -  ki  V2  +  k2  V5  (31) 


where : 


•  k^  is  a  matrix  of  constants  containing  the  number  of  rows  equal 

to  the  number  of  class-Y  voltage  sources  and  the  number  of  columns 
equal  to  the  number  of  non-scannable  class-2  elements. 

•  Y.^  is  a  matrix  of  constants  containing  the  number  of  rows  equal 

to  the  number  of  class-Y  voltage  sources  and  the  number  of  columns 
equal  to  the  number  of  non-scannable  class-5  elements. 
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•  EY  Is  a  vector  composed  of  class-Y  voltage  sources. 


Also,  since  any  resistor-current-dependent  current  source  must  depend  on  a 
resistor  tree  branch  current  or  a  resistor  link  current,  there  arises 


JX  -  k.  I0  +  k.  I. 

3  2  4  5 


where 


•  kj  is  a  matrix  of  constants  containing  the  number  of  rows  equal 
to  the  number  of  class-X  current  sources  and  the  number  of 
columns  equal  to  the  number  of  non-scannable  class-2  elements 

•  is  a  matrix  of  constants  containing  the  number  of  rows  equal 

to  the  number  of  class-X  current  sources  and  the  number  of  columns 
equal  to  the  number  of  non-scannable  class-5  elements 

•  JX  is  a  vector  composed  of  class-X  current  sources 

If  Equations  (31)  and  (32) ,  together  with  V?  ■  R__I_  ant*  "  ^SSV5  are 
substituted  into  Equations  (29)  and  (30) 


‘22*2  +  B25V5  + 

B2y  |klB22I2  +  k2V5|  * 

-B24V4 

'55V5  *  B25  T2 

-bx5t  k3i2  +  V55v5| 

-  B35 

The  last  two  equations  may  be  consolidated  as 

f  (R22  *  VlR22>  <B25  +  B2yk2>  1  f  *2! 


("B25  "  Bx5  k3}  (G55  '  Bx5  k4G55 ^  1  LV 


-B24V4  -  B27E7 


T  T 

B35  13+B85  J8 


If  the  large  matrix  on  the  left  side  is  called  MRG,  then 
m  r  -B2«V4  -  B27E7  1 


T  T 

B35  I3  +  B85  J8 


and  the  resistive  quantities  and  can  be  determined  without  computational 
delay.  Note  that  since  all  four  "k"  matrices  are  constrained  to  be  constant. 
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Che  linearly  dependent  source  feature  Itself  will  not  require  any  more  than 
one  inversion  of  MRG.  If  any  variable  resistors  are  present  in  a  network, 
MRG  must,  of  course,  be  inverted  at  each  time  step.  In  addition,  the  extra 
row  and  column  that  was  added  to  the  B  matrix  of  figure  2  by  the  linearly 
dependent  sources  will  add  terms  to  the  equations  used  in  the  solution  of 
capacitive,  inductive  and  some  source  quantities.  Specifically,  Equations 
(5),  (3),  (4)  and  (8)  become 

X4  "  BUT  h  +  B24T  ;2  +  B34T  *3  +  "Si"  J8  +  Bx4T  ^  (5,) 


V,  -  -B-.V.  -  B,,  V,  -  B,,  V,  -  B„  E,  -  B,  EY 

3  34  4  35  5  36  6  37  7  3y 


V8  '  -B84V4  -  B85  V5  '  B86  V6  '  B87  E7  '  V  " 

X7  ■  BX7T  X1  +  B27T  X2  +  B37T  X3  +  B87T  J8  +  Bx7T  JX 


and  the  new  relations 


™  -  -Bx4V4  -  Bx5V5  -  Bx7E7  *  BxyEY  (35> 

tY  ’  +B2yT  h  +  B3yT  X3  +  B8yT  J8  +  BxyT  JX  <36> 

now  exist. 

A  restriction  must  be  placed  on  these  sources  based  on  the  content  of 
section  2.5.  The  B  matrix  of  figure  5  transforms  Equations  (9)  and  (10)  into: 


-B..  V.  -  R.,  E,  -  B-  EY 
14  4  17  7  ly 


X6  '  B36T  X3  +  B86T  V  Bx6T  3X 


(10') 


Now,  additional  time  derivatives  EY  and  JX  are  required  whenever  B..  or 
Bx6  *  °-  y 


.  Differentiation  of  Equations  (31)  and  (32)  would  involve  quantities  1^, 
V^,  I^,  and  t^.  The  formulation  contains  no  provisions  for  these  quantities, 
and  there  is  no  way  the  user  could  know  them  to  supply  them  as  input  data. 
SCEPTRE  will  automatically  check  for  the  existence  of  non-zero  B^  or  B  ,  and 
terminate  the  run  with  a  diagnostic  message  when  they  occur.  Theysituation 
can  be  recognized  by  the  presence  of  any  circuit  loop  composed  solely  of 
capacitors  and  at  least  one  linearly  dependent  voltage  source,  or  the  presence 
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of  any  circuit  cut  set  composed  solely  of  inductors  and  at  least  one  linearly 
dependent  current  source. 

2.3.5  INTEGRATION  ROUTINES 

Four  integration  routines  are  optionally  available  for  use  in  SCEPTRE. 
Two  of  these,  RUK  and  TRAP,  were  available  in  PREDICT  and  have  been  only 
slightly  modified.  The  third  routine,  called  XPO,  was  developed  at  the  IBM 
Scientific  Center  in  Palo  Alto,  California,  by  Dr.  R.  Warten  and  Mr.  M.  Fowler 
and  adapted  for  use  in  SCEPTRE.  Studies  to  date  indicate  that  XPO  is  usually, 
although  not  always,  faster  than  the  other  methods.  For  that  reason,  this 
routine  will  always  be  used  unless  the  user  explicitly  requests  otherwise 
in  the  RUN  CONTROL  section  of  any  run.  The  fourth  integration  routine, 
IMPLICIT,  is  discussed  in  subsection  2. 3. 5. 5.  This  routine  has  advantages 
when  the  spread  between  the  largest  and  smallest  time  constants  in  a  network 
to  be  analyzed  is  very  large. 

2. 3. 5.1  RUK  Formulas  and  Variable  Step  Size  Control 

The  well-known  Runge-Kutta  fourth-order-accuracy  formulas  (Ref.  2) 
for  the  numerical  solution  of  the  differential  equation 

y  -  f  (t,  y) 

are  given  by 

y  (t  +  h)  -  y  (t)  +  1/6  k^  +  2k2  +  2k3  +  k^  (37) 

where 

k3  *  h  f  t,  y  (t) 

k2  -  h  f  J  t  +  y  (t)  +  f  kx 

k3  "  h  f  \l  +  2’  y  (t)  +  2  k2  | 

k^  -  h  f  t  +  h,  y  (t)  +  kj 

Variable  step  size  control  (Ref.  3)  is  achieved  by  computing 
k5  -  h  f  t  +  h,  y  (t  +  h)j 
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and 


E  -  k,  -  2k,  -  2k.  +  3k. 
13  4  5 


(38) 


The  above  formulas  easily  generalize  Co  systems  of  differential  equations.  In 
which  case  y  and  E  become 

y  (t)  -  [yx  (t> ,  — ,  yn  (t)] 

E  "  <ei’  .  en> 

Now  let 

u.  i  u2»  »  3. 2  £  0 

Set 

Uk  ‘  “,*“2  1  \  (t  +  h)  I 

\  -  ll  +  h  I  yk  (t  +  h>  I 

If  |  |  >1.5  for  some  k,  the  Integration  step  h  is  halved,  the  in 

dependent  variable  Is  restored  from  t  +  h  to  t,  and  the  values  of  y  and  y  are 
restored  to  the  values  at  time  t. 


If  |  |  >  0. 75  for  some  k,  and  |  e^  |sl.5  for  all  k,  the  current 

Integration  step  is  accepted,  but  the  step  size  is  halved  for  succeeding  steps. 


If  |  e^  1^0.75  for  all  k, 


the  step  size  is  unaltered. 


If  |  ejJ  s  for  all  k,  and  |  e^|<L^  for  at  least  one  k,  a  doubling  in¬ 
dicator  is  activated.  Actual  doubling  isKdelayed  for  seven  time  steps. 
Halving  always  takes  precedence  over  doubling;  thus,  anytime  a  halving  signal 
is  received,  the  step  size  is  halved  and  doubling  is  delayed  for  at  least 
seven  steps.  Similarly,  after  successful  doubling,  another  seven  steps  must 
elapse  before  the  step  size  can  be  doubled  again. 


Recommended  choices  for  u^,  u2>  1^,  and  lj  are  as  follows.  If  absolute 
error  control  is  desired, 

set 

u^  *  0.0075  U2  ■  0 

ll  -  0.00005  i2  -  ° 
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If  relative  error  control  Is  desired, 
set 

u1  -  0.005  u2  -  0.005 

lx  -  0.00005  12  -  0.00005. 

If  smaller  step  sizes  are  desired  than  the  ones  yielded  by  the  above  set¬ 
tings,  the  0.0075  and  0.00005  settings  should  be  reduced  by  a  factor  of  32. 
This  will  yield  half  the  previous  step  sizes. 


2. 3. 5. 2  Modified  Trapezoidal  Integration  (TRAP) 
2. 3. 5. 2.1  Method 

Given  the  differential  equation 


y(t) 


f  (t,  y  (t)),  y  (0) 


consider  the  following  numerical  integration  scheme, 

y  (t  +  |  h)  “  y  (t)  +  j  y  (t) 


yp  <c  +  f  h> 


t  +  f  h»  yp  (t  +  f  h) 


t  +  j  h,  y  (t)  +  |  y  (t) 


yc  (t  +  h)  -  y  (t)  +  j 


2  y  (t)  +  y„  (t+7  h) 


yc  (t  +  h) 


t  +  h,  yc  (t  +  h) 


Step-size  control  is  achieved  by  computing 

E  "  §  h  I  yp  (t  +  \  h)  -  y  (t)  | 


The  magnitude  of  E  inuicates  any  changes  to  be  made  in  step  size. 


(39) 

(40) 

(41) 

(42) 

(43) 


2. 3. 5. 2. 2  Truncation  Error 


Inasmuch  as  y  ■ 

ft  +  y 

yT  (t  +  h)  - 

y  (t) 

f  and  the  true  solution 

y  2  .. 

+  hy  (t)  +  |  y  (t) . 


yT 


is  approximated  by 
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one  obtains, 


yT  (t  +  h)  -  y  (t)  +  h  y  (t)  +  (ffc  +  y  fy> 


(44) 


On  the  other  hand 
f 


t  +  y  (t)  + 1  y  (t) 


whence  Equation  (41  )  becomes 

yc  (t  +  h)  1  y  (t)  +  y 


f  (t,  y  (t))+  ^  ffc  +  &  fy 


2y  (t)  +  y  (t)  +  ^  ft  +  |  y  f 


Subtracting  Equation  (41  )  from  (44)  yields 


y_  (t  +  h)  -  y  (t  +  h)  A 
i  c 


(4-0  * 


•  r  »»  •  c 

y  fy  ■  r  y  fy 


which,  in  turn,  yields  the  approximate  local  truncation  error. 


(45) 


From  Equations  (40)  and  (45)  one  gets 

•  f  .  3h  .  ■  f  \  .  3h  r*  .  h  •  r- 

yp  (t  +  y  (t)  -  -J  ft  2yfy 


whence 


2h 

3 


y  (t  +^)  -  y  (t) 


9  h2 

4  h  ft+T-yfy 


(46) 


One  notes  that  if  ffc  *  0,  then  Equation  (46)  is  a  good  representation 

of  the  local  truncation  error,  neglecting  higher  order  terms. 

Experience  indicates  that,  for  systems  of  equations  of  the  form  Y  = 

A  (y)  Y  +  B  where  A  is  piecewise  constant,  the  method  as  well  as  step-size 
control  is  adequate. 


2. 3. 5. 2. 3  Stability 

Consider  the  differential  equation 
y  *  -ay.  y(0)  -  yQ,  a  >  0 


The  numerical  integration  scheme  given  by  Equations  (39)  through  (42) 
yeilds 

yc  (nh)  -  1  1  -  ah  +  --jP-  ■  j  y  (0)  (47) 
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This  Is  easily  established  by  induction.  Thus,  for  numerical  stability 

C  h)  2 

it  is  necessary  and  sufficient  that  |  1  -  ah  +  —  |<1.  This  yields  0  < 

ah  <  6  as  shown  in  figure  6. 

Inasmuch  as  |  ah  |  <  6  is  necessary  and  sufficient  for  numerical  stability, 
we  say  that  the  modified  trapezoidal  method  has  a  stability  radius,  r,  equal 
to  6.  Moreover,  from  Equations  (39)  and  (40)  we  see  that  the  method  requires 
two  derivative  evaluations  (passes)  per  integration  step.  Thus,  the  pass 
number  p  associated  with  the  method  is  2. 

The  ratio  r/p  is  a  measure  of  a  method's  efficiency  in  the  sense  of  minimum 
number  of  integration  steps  per  given  time  interval. 


ah 


Figure  6.  Root  Locus  of  TRAP  Characteristic  Equation 


The  higher  the  r/p,  the  fewer  steps  are  required  for  a  given  a. 
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The  following  list  shows  r,  p  and  r/p  for  a  few  representative  methods: 


METHOD  r 

Modified  TRAP  6 

Euler  (Ref.  4)  2 

Trapezoidal  (Ref.  5)  2 

NIDE  (1PASS)  (Ref.  6)  0.8 

Runge-Kutta  (Fourth-order)  2.78 

NIDE  (2PASS)  1 

Hamming  (Ref.  7)  0.85 


£ 

2 

1 

2 

1 

4 

2 

2 


l/£ 

3 

2 

1 

0.8 

0.7 

0.5 

0.425 


2. 3. 5. 2. 4  Step-Size  Control 

As  mentioned  previously,  the  quantity 

e2  "  r  y p  <*  +  r>  -  y  (t) 

is  used  for  step-size  control. 

(1)  The  method  is  essentially  the  same  as  for  RUK.  Let  u^,  u^,  1^,  12 
be  real  numbers  2  0. 

Compute  U  -  u..  +  y£  (t  +  h)  ,  L  -  1^  +  y  ^  (t  +  h) 

If  Ej  >  0.75  U,  the  step  size  is  reduced. 

If  E2  <  L,  the  step  size  is  increased. 

If  L  S  Ej  £  0.75  U,  the  step  size  is  not  altered. 

(2)  Choice  of  u^,  u2>  1^  12- 

Recall  that  E2  attempts  to  represent  the  local  truncation  error.  For 
absolute  error  control  set 

u^  ■  0.01  ,  u2  ■  0 

lx  -  0.0005  ,  12  "  0 
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For  relative  error  control  set 


■  0.001 
0.00005 


u2  -  0.01 

0.0005 


‘1  -  .  ’  2 

The  above  values  work  well  in  practice. 

Reducing  the  above  setting,,  by  a  factor  of  4  will  reduce  the  step  sire 
by  a  factor  of  2. 

2. 3. 5. 3  Exponential  Integration — (XP.Ql 

2. 3. 5. 3.1  Method 

Given  the  differential  equation 

y  (t)  -  f  1 1 ,  y  (t)  | ,  y  (0)  -  y0  * 

this  integration  routine  gives  the  computed  solution 


Xh  _  a 

y  (t  +  h)  -  y  (t)  +  yA  (t>  h  +  h  *P  (t) 

except  for  two  cases  that  will  be  covered  later.  In  equation  (48) 


yA  (t)  - 


y  (t)  -  y  (t  -  hQ) 


(48) 


»here  t  -  hQ  is  the  last  point  computed; 


yp  (t) 


y  (t)  -  yA  (t); 


y  (t)  /  yp  (t) 


(49) 

(50) 


Since  y  <t)  cannot  be  computed  explicitly  by  the  program,  we  take  a  small 
Euler  integration  step  and  approximate  y  (0  y 

y  (t  +  6  )  ■  y  (t)  +6  y  (t) 

6 


y  (t  +  6)  -  f  I  t  +  6,  y6  <t  +  5)  |  ^ 


V°  - 


y6  (t  +  6  )  -  y  CO 


,  0  <  6  *  4 


a  „  4f  v  and  y  (t)  have  the  same  sign, 

rhe  first  exception  occurs  if  y ^  tw  anQ  yp  v 
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Then  sec 


yA  <0  -  0, 


which  results  In  Equations  (49)  and  (SO)  becoming 

yp  (t)  -  y  (t)  (52) 

X  -  yfi  (t)  /  y  (t). 

This  provides  a  better  approximation  to  the  solution  of  the  equation 
y  ■  X  y  +  b  without  decreasing  the  overall  effectiveness  of  the  method. 

The  second  exception  occurs  If  X  ,  computed  either  by  Equations  (50)  or  (52) 
is  non-negative.  Then  the  term 


in  Equation  (48)  is  replaced  by  1  +  — j.  thus  changing  Equation  (48)  to 

y  (t  +  h)  -  y  (t)  +  yA  (t)  h  +  (1  +  ^— )  hyp  (t)  (53) 

It  can  be  shown  with  some  effort  that  Equation  (53)  Is  equivalent  to  a  trap¬ 
ezoidal  type  Integration  method. 


2. 3. 5. 3. 2  Truncation  Error 

x  -i 

6  “1 

Noting  that  — - - ,  when  expanded  in  a  Taylor  series,  becomes 

.  l+i  +  i!  +  0  (x3, 

Equation  (48)  can  be  written 

h2 

y  (t  +  h)  -  y  (t)  +  yA  (t)  h  +  yp  (t)  h  +  -j  Xyp  (t)  + 

1,3  t 2  r„  <-t>  +  0  (h4>  (54) 

6 

2  3 

-  y  (t)  +  y  (t)  h  +  J-  y&  (t)  +  |-Xy6  (t)  +  0  (h4) 
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Referring  to  Equation  (51) ,  note  that 


y*  (t)  -  y  (t)  +  f 


y  (t)  -  y  (t)  ^  +o  (62) 


(55) 


Substituting  Equation  (55)  into  equation  (54)  yields 

2  2 

y  (t  +  h)  -  y  (t)  +  y  (t)  h  +  y  (t)  j-  +  £-6 


y‘(t)  - 


y  <o  |S 


+  o  (hJ) 
56) 


It  is  noted  that  the  use  of  either  exception  to  the  principal  Equation  (48) 
still  yields  the  truncated  expression,  Equation  (56). 

_  Since  6  5  the  method  is  second  order  exact;  i.e.,  the  error  is  of  order 
h  .  This  characteristic  will  be  used  in  Equation  (57)  in  the  next  section. 


2. 3. 5. 3. 3  Step-Size  Control 


For  t  £  t  +  ££t  +  h,  let  y^,  (t  +  £  )  be  the  true  solution  to  the 
differential  equation,  let  y  (t  +|)  be  the  computed  solution,  let  yg  (t  +§  ) 
be  an  expression  obtained  bycdifferentiation  of  yc  (t  +5)  with  respect  to 

'  yA  (t)  +  (1  +  xS)  y  (t),  X2o 


ye  (t  +  £  )  * 


( yA  (t)  +  e  \%y  (t) 


X<0 


and  let  y  (t  +  £)  be  obtained  by  substituting  y  (t  +£)  into  the  differential 
equation. 


Ideally,  step  control  should  be  based  on  the  expression 
-h 


■ru 

J  0 


yT  (t  +  % )  -  yg  (t  -  O 


d  % 


however,  y^  (t  +  %  )  is  not  available  except  at  §  ■  0. 

Making  use  of  the  fact  that 


yc  (t+5) 


f  1 1  +  5  ,  yT  (t  +  %  )  +  |  yc  (t  +  §  )  -  yT  (t  ♦  5  )  |  j 


yT  (t  +§)  +  0  (5J) 


(57) 


we  obtain 
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h  h 

E0  -  j  [*c  <t+5> -*.<'+«>]  +y  [yt  (t  +5)  - 

J  0  0 

yc  (t+  S>]  d  5 

h 

^  (t  +?)  -  y.  (t  +5)  d  5  +  0  (h4) 

c  e  j 

0 

•  •  »* 

If  yc  (t  +5  )  -  yg  (t  +  §)  does  not  change  sign  for  0  s  ^  s  h, 

|  EQ|  i  h  |yT  (t  +  h)  -  yg  (t  +  h)  |  s  h|  yc  (t  +  h)  -  yg  (t  +  h)  +  0  (h4) 

Thus,  for  small  h, 

E  ■  h  |  yc  (t  +  h)  -  ye  (t  +  h)| 

yields  an  estimate  of  the  truncation  error. 

Let  u^,  u2>  1^,  12  be  real  numbers  2  0.  Compute 

U  -  ^  +  u2  |  y£  (t  +  h)  J 

L  -  lx  +  12  |  yc  (t  +  h)  I 

If  |  E  |>U,  the  step  size  is  reduced.  If  |  Ej  <  L,  the  step  size  is  increased. 

If  L  s  E  £  U,  the  step  size  is  not  altered. 

For  absolute  error  control,  set 

u^  ■  0.0075  ,  u2  *  0 

1L  -  0.0002  ,  i2  -  ° 

For  relative  error  control,  set 

u^  ■  0.005  ,  u2  ■  0.005 

1  -  0.0001  ,  l2  -  0.0002 
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2. 3. 5. 4  General  Absolute  and  Relative  Error  Criteria 


All  three  explicit  Integration  methods  described  In  this  section  are 
approximate  methods  that  are  used  to  integrate  differential  equations  and, 
thereby,  update  the  state  variables  at  each  time  step.  Each  of  these  methods 
has  an  automatic  control  that  allows  step  size  to  be  increased  and  decreased 
during  the  transient  problem.  The  method  of  control  used  compares  some  func¬ 
tion  of  the  step  size  and  the  derivatives  of  the  state  variables  to  a  quantity 
that  serves  as  a  standard  as,  for  example,  in  Equation  (43).  The  standard 
may  be  a  constant  or  a  variable  function  of  the  state  variable.  If  the  former 
is  used,  it  is  termed  an  absolute  error  criterion;  if  the  latter  is  used,  it 
is  termed  a  relative  error  criterion. 

Consider  the  situation  in  which  two  state  variables  are  very  different  in 
size.  Let  y^  (t)  ■  1  and  (t)  ■  100.  For  all  practical  purposes,  it 
is  usually  unnecessary  to  integrate  the  derivatives  of  these  state  variables 
to  the  same  accuracy.  If  an  absolute  error  criterion  is  used,  this  is  just 
what  is  done,  and  the  step  size  may  be  unnecessarily  inhibited.  If,  however, 
a  relative  error  criterion  is  used,  effectively  a  larger  error  will  be  tol¬ 
erated  for  the  integration  of  y.  (t)  and  a  larger  step  size  will  be  permitted. 

In  general,  then,  it  would  be  best  for  the  user  to  use  relative  error  control 
when  large  values  of  state  variables  (capacitor  voltages  and  inductor  currents) 
are  expected. 

Relative  error  control  is  programmed  with  all  three  explicit  methods, 
but  the  user  may  easily  modify  the  relative  controls  or  enter  absolute  con¬ 
trols  (subsection  2.2.10  of  Volume  I).  If  larger  numbers  are  used  tor  the  u^  or 
u2  entries,  the  solution  process  will  be  less  likely  to  halve  any  particular 
solution  step  size;  smaller  numbers  would  increase  the  likelihood  of  reduced 
step  sizes.  If  larger  numbers  are  used  for  the  1^  and  1^  entries,  the  solu¬ 
tion  process  is  more  likely  to  increase  any  particular  solution  step  size; 
smaller  numbers  make  Increased  step  sizes  less  likely.  The  user  should 
realize  that  any  increase  in  solution  speed  that  may  result  from  adjustment 
of  these  numbers  must  necessarily  come  at  the  expense  of  integration  accuracy. 

2.3.5. 5  Implicit  Method 

The  integration  methods  described  in  subsections  2. 3. 5.1,  2. 3. 5. 2 
and  2. 3. 5. 3  are  all  explicit  in  form  and  can  be  used  interchangeably  within 
the  mathematical  formulation  of  SCEPTRE  without  difficulty.  If,  however,  any 
implicit  method  is  to  be  used,  whether  single  step  or  multistep,  an  additional 
computational  step  is  necessary.  This  step  is  discussed  in  the  following 
paragraphs  and  in  Appendix  A  . 
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•■wv .  t.yj;  „ 


2. 3. 5. 5.1  Basic  Implicit  Format 

A  generalized  form  of  all  implicit  integration  methods  can  be 

written  as 

P  P 

Y (N+l)  -  £  a.  Y(n-i)  +  h  £  b.  Y(n-i)  (58) 

i-0  i— 1 


where  Y  is  the  vector  of  system  state  variables,  Y  is  the  vector  of  state 
variable  derivatives,  n  is  the  step  number,  h  is  the  step  size  and  the  a^, 
b ;  are  suitably  chosen  constants.  Multistep  methods  are  introduced  if 
P>0«  The  simplest  derivative  form  of  Equation  (58)  is 

• 

Y (n+l)  -  Y(n)  +  h  Y(n+1) 


which  is  commonly  referred  to  as  the  implicit  or  backward  Euler  technique. 
If  an  mth  order  system  of  differential  equations  is  to  be  solved  by  this 
method,  the  following  generalized  matrix  equation  will  result: 


r  aYi 

8Y1 

1  “ hg  aYx 

(Yl,...Ym,t)  . 

.  •  -hg  gY 

•  m 

(Yl,...Ym,t) 

•  • 

•  9Y 

•  ♦ 

* 

'h8_5*7 

(Yl,...Ym,t)  . 

.  .  1-hg 

°  m 

-  (Yl...Ym,t) 

"  -F^  (Yl,..Ym,t) 

~F  (Yl,..Ym,t) 
m 


(59) 


The  iteration  implied  in  Equation  (59)  is  carried  out  to  convergence  at  each 
time  step.  The  k  superscripts  here  indicate  the  kth  approximation  to  the 
final  value  at  convergence  at  each  step,  g  is  a  constant  that  depends  on  the 
order  of  integration,  and  F^  is  a  function  of  the  ith  differential  equation, 
the  step  size  and  past  values  of  the  mth  state  variable.  Sparse  matrix 
techniques  will  be  applied  to  the  operation  implied  by  Equation  (59)  when 
large  problems  are  encountered. 
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2. 3. 5. 5. 2  The  Jacobian 

The  Y^  terms  in  Equation  (59)  are  readily  available  from  the  basic 
program  formulation;  but  the  general  partial  derivative  term 


l£i£m,  l£j&m 


is  not.  To  determine  what  is  really  needed,  it  is  desirable  to  frame  the 
entire  derivation  of  the  symbolic  Jacobian  in  terms  of  the  general  matrix 
equation 

Y  -  AY  +  BU  +  NU  (60) 


Since  the  desired  quantities  are  the  general 


it  can  be  seen  from  partial  differentiation  of  Equation  (60)  that  these  are 
contained  in  the  matrix  A.  Hence  the  convenience  of  the  general  notation  is 
given  in  Equation  (60). 


What  now  follows  is  the  construction  in  symbolic  form  of  the  general  matrix  A 
in  terms  of  the  mathematical  formulation  that  was  derived  in  subsection  2.3.1. 
Begin  with 


C44V4 


’14 


V  B24 


T  T 

*2  +  B34 


*3  +  B84 


(5) 


and 


V3  “  _B34  V4  "  B35  V5  "  B36  V6  “  B37  E7 


(4) 


Substitute  Equation  (9)  and  (13)  into  equation  (5)  for 

C44  \  -  BU  T  CU(-B14VB17V  +  B24  T  V1  <-B2«VB27E7 
'B25  R55  B35  1 3  '  B25  R55  B85  J8*  +  B34  *3  +  B84  J8 


(61) 
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Use  Equations  (24),  (10),  (23)  and  (15)  into  Equation  (4)  for 


L33B3  +  l36  <B36  "  *3  +  B86  "  V 


"  "B34V4  "  B35  0 


-1 


-B25  T°22  (B24V4  +  B27E7>  +  B35  T  *3  +  B85  ' 


-  8 36  a63  +  L66  B36T)  *3  '  B36  L66  B86  '  G8  ‘  B37E7 


Equation  (61)  can  be  manipulated  to  yield 

<C44  +  B14  1  CnB14>  E4  '  -B14  "  CUB37  E7  *  B24  T  ‘‘  (B24  V4 


+  B27E7  +  B25R55B35  T  *3  +  B25R5SB85  '  J8, 


T  T 

+B34  I3+B84  J8 


Equation  (62)  yields 
[L33 

"o  "  (- 

)  -B37E7  - 


+  L36  B36  T  +  B36  L63  +  B36  L66  B36  "]  *3 


'  -B34V4-B35 


B25  G22B24V4  "  B25  G22B27E7  +  B35  *3  + 


B85  T  J8 


T  •  x  ’ 

L36B86  J8  "  B36L66B86  J8 


If  only  the  coefficients  of  V4,  I3,  V4  and  I3  (the  stat^  variably 
state  variable  derivatives)  are  retained,  Equations  (61)  and  (62) 
considerably  simplified  and  can  be  written  in  matrix  notation  as 


-1 


■  “ 

V. 

M  0 

4 

c 

• 

*3 

1 - 

0 

L  ... - 

(62) 


(61) 


(62) 


s  and  the 
become 
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<-B24  T  ”r  <B34  T  -  B24  *  “r  '  B25R55B35 

(B35M0  B25  G22B24  '  B34)  ('B35MG  B35  ’ 


where  Mc  -  +  >u  Cu 

and  -  L33  +  L36B36  T  +  B36L63  +  BJ6L66B36  T 

The  matrix  appearing  on  the  right  side  of  Equation  (63)  is  now  in  a  form 
corresponding  to  the  first  term  on  the  right  side  of  Equation  (60).  This 
matrix  is  the  symbolic  form  of  the  general  A  matrix  that  is  needed  to  use 
implicit  integration  (Ref.  8,  9)  with  the  basic  program  formulation. 

Another  method  is  available  to  construct  the  Jacobian  that  is  based  on  a 
numerical  rather  than  a  symbolic  approach.  Each  time  that  the  Jacobian  is 
to  be  re-evaluated,  m  calls  are  made  to  SIMUL  8  to  compute 


F(Y1,  TJt  ...I  .Y  ♦«  Y  ,  •••  Ym,  C) 


Approximations  to  the  desired  partials  are  than  obtained  by 


9Yi 


;<j> .  * 


2.4  DC  SOLUTIONS 

SCEPTRE  offers  five  DC  solution  modes,  Initial  Conditions,  Sensitivity, 
Monte  Carlo,  Worst-Case  and  Optimization.  Any  one  of  the  five  may  be  used 
alone,  or  as  a  source  of  initial  conditions  for  a  Transient  or  AC  run. 
Alternatively,  the  user  may  supply  initial  conditions  data  himself  as  entries 
to  a  transient  or  AC  run  without  the  help  of  the  DC  options.  One  of  the  DC 
options,  called  Initial  Conditions,  takes  part  in  all  DC  solutions.  It  may 
be  called  separately,  and  if  any  other  DC  option  is  called,  that  option  uses 
two  or  more  passes  of  the  Initial  Conditions  calculation  to  produce  its 
results.  The  Initial  Conditions  solution  (and  thus  the  other  DC  options 
also)  may  use  either  the  Newton-Raphson  or  implicit  method.  (See  App.  A,  Vol.  I) 

2.4.1  INITIAL  CONDITIONS  SOLUTION  BY  NEWTON-RAPHSON  METHOD 

2. 4. 1.1  Technique  Description 


Many  practical  circuits  require  computer  solution  of  the  initial 
conditions  prevailing  at  the  start  of  the  transient  (time  before  the 
transient  solution  can  begin.  These  values  can  always  be  determined  by  a 
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separate  transient  run  In  which  all  forcing  functions  are  held  at  the  values 
for  time  -  t0.  However,  an  alternate  procedure  based  on  an  Iteration  technique 
using  independent  variables  other  than  time  was  considered  desirable.  This 
procedure  presents  the  advantage  of  economy  of  machine  time  on  all  circuits 
for  which  convergence  occurs.  This  section  will  describe  the  formulation  of 
this  portion  of  SCEPTRE,  which  is  completely  independent  of  the  transient 
formulation. 

A  similar  derivation  for  the  adjoint  calculation  is  given  in 
Appendix  P. 

If  an  L  tree  is  set  up  based  on  the  preference  order  E,  L.  R,  C  a 
general  B  matrix  may  be  set  up  according  to  the  procedure  outlined  in  sub¬ 
section  2.2.  The  zero-valued  submatrices  arise  from  the  L  tree  and  the 
preference  order*.  The  resulting  B  matrix  is  shown  in  Table  II. 


TABLE  II 

GENERAL  B  MATRIX  FROM  THE  L  TREE 


4 

5 

6 

1 

7 

1 

B14 

B15 

B16 

B17 

2 

0 

B25 

B26 

B27 

G3 

0 

0 

B36 

B37 

8 

B84 

B85 

B86 

B87 

9 

B94 

B95 

B96 

B97 

1 

1 

B04 

B05 

B06 

B07 

The  following  equations  (among  others)  arise  from  this  B  matrix  if 
vectors  V,,  I.,  1^  and  submatrix  B^  are  assumed  to  be  zero.  These  assump¬ 
tions  are  based  on  the  known  final  values  of  V^,  1^,  and  1^  for  the  initial 
condition  problem  and  the  absence  of  any  current-source  capacitor  cut  sets 
(see  the  restrictions  in  subsection  2, 4. 1.4). 


For  example,  the  submatrix  B24  must  always  be  zero  since  non-zero  entries 
in  it  could  only  arise  when  resistor  links  close  loops  containing  capacitor 
tree  branches.  The  preference  order  prohibits  this  possibility. 
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V9  +  B95  V5  +  B97E7  -  0  (65) 

V2  +  B25  V5  +  B27  E7  "  0  <66> 

h  -  B25  T  h  -  B85  T  J8  -  B95  '  J9  '  B05  "  J0  *  0  <67'> 


To  get  all  the  variables  of  Equation  (67)  in  terms  of  Vg,  V^,  V^,  and 
effective  sources,  the  following  substitutions  are  made: 


*5  ‘  G55  V5 


*2  "  C22  V5 


J9  “  G99  V9  +  ^9 
where: 

•  a  is  a  matrix  containing  the  number  of  rows  equal  to  the  number 
of  secondary  current  sources  and  the  number  of  columns  equal  to 
the  number  of  primary  current  sources 


•  Gj.^  and  G ^  are  diagonal  matrices  containing  only  conductances 

#  G99  is  a  diagonal  matrix  containing  only  diode  and  transistor 
junction  conductances 


•  Q  terms  are  described  in  Appendix  II 
Then,  Equation  (67)  becomes 


C55  V5 


25 


T  G22  V2  '  B85  T  J8  - 
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Equations  (65),  (66),  and  (67)  may  be  designated  as  (V9,  V2,  V5) , 
f2  (v9*  V2*  v5^  and  F3  (v9*  v2*  v5^  »  respectively.  If  the  basic  Newton- 
Raphson  method  (see  Appendix  D  )  is  applied  to  Equations  (65) ,  (66) ,  and 
(67): 


F1  (V9*  V2*  V  + 

^1 

/  \;  \/  17  \ 

\\7  4- 

dV9 

W91  ^  2* 

AV9  + 

dFl 

—2  <V  V  V 

AV2 

+ 

dFl 

—5  <V  V2*  V 

av5  -  0 

F2  (V  V2»  v5>  + 

*F2 

SV9 

(V9,  v2,  v5) 

av9  + 

3F2 

-Sv;  <V  V  V 

AV2 

■f 

*F2 

—2  <V  V  V 

av5-o 

F3  <V  V2 ’  V5)  + 

3F3 

<v9,  v2,  V,) 

\tr  4. 

9V9 

^V9  + 

BF3 

“SVJ  (V  V  V 

^V2 

+ 

dF3 

8  V,  <V  V  V 

^5 

-  0 

or 


F1  (V  V  V 


F2  <v9,  v2,  v5) 


f3  (V9,  V2,  V5) 


+  z 


av„ 


avp 


(68) 
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where  the  Jacobian  is 


aFl 

(V 

*F1 

*Fi 

*V9 

>  •  ^  t 

v2*  V5) 

aV2 

vv9 1 

V  V 

oy5 

9F2 

dF2 

/IT 

it  17  \ 

&F2 

SV9 

(V 

V2’  V 

SV2 

V  | 

V2*  V5 

5V5 

dF3 

(V9, 

V  V 

aF3 

/v 

V  v  \ 

SF3 

dV9 

dV2 

v*g  t 

V  v5) 

bv5 

(VQ,  V,,  V.) 


(VQ,  V.,  V.) 


(VQ,  V0,  V,) 


(69) 


If 


•^V9  -  V9  (n  +  1)  -  V9  (n) 
AV2  -  V2  <n  +  1)  -  V2  (n) 
J^V5  -  V5  (n  +  1)  -  V5  (n) 


where  n  is  used  to  designate  the  results  of  the  n£b  iteration  pass,  then 
Equation  (68)  becomes 


—  “ 

F1  (V9 ’  V2*  V 

V9  (n  +  1)  -  V9  (n) 

f2  (v9,  v2,  v5) 

+  Z 

V2  (n  +  1)  -  V2  (n) 

f3  (v9,  v2,  v5) 

V5  (n  +  1)  -  V5  (n) 

leading  directly  to 


V9  (n  +  1) 

V9  (n) 

Fx  (V9,  V2,  V5) 

V2  (n  +  1) 

- 

/-S 

c 

CM 

> 

-z-1 

F2  <V  v2  *  v5> 

V5  (n  +  1) 

V5  (n) 

f3  (v9,  v2,  v5) 

_ 

_ 

—  — 

(71) 
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Equation  (69)  may  be  written  more  explicitly  if  the  indicated  differentiations 
are  performed  on  Equations  (65) ,  (66),  and  (67)  to  obtain: 


so  that  Equation  (71)  becomes 

V9  (n  +  1)  V9  (n) 
V2  (n  +  1)  -  V2  (n) 
V5  (n  +  1)  V5  (n) 


<v9,  v2, 
(V9,  v2, 


(72) 


Equation  (72)  may  be  written  in  more  convenient  form  as  follows  (see  also 
Appendix  E )  : 


Equation  (72')  signifies  a  computational  sequence  as  follows: 

Quantities  Vg,  V2,  and  V5,  the  vectors  of  voltages  across  primary  current 
sources,  resistor  link,  and  resistor  tree  branches  respectively,  each  have 
some  assumed  value  (usually  zero)  to  begin  the  computation  at  n  *  0.  All 
members  of  the  right  side  of  Equation  (72')  that  may  be  voltage  dependent 
are  updated.  The  left  side  of  Equation  (72')  is  then  computed  and  the  first 
iteration  (n  =  1)  is  complete.  The  right  side  of  Equation  (72')  is  re-evaluated 
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on  the  basis  of  the  results  of  the  first  iteration  and  the  left  side  is 
again  computed,  thereby  completing  the  second  iteration  (n  ■  2).  This 
process  is  repeated  up  to  100  times.  After  any  iteration,  if 

l  v<„  + 1)  -v(„)li0-001  l  v<„  ♦  i)  l +  °-0001  <73) 


is  satisfied  for  the  set  of  all  voltages  in  Equation  (72) ,  convergence  is 
considered  to  have  occurred  and  the  procedure  is  terminated.  If  after  100 
iterations  Equation  (73  )  is  not  satisfied,  a  diagnostic  message  will  be 
printed  indicating  that  convergence  has  not  occurred.  Experience  to  date 
has  shown  that  the  Newton-Raphson  procedure  will  converge  for  most  circuits 
in  less  than  30  iterations. 


Convergence  of  the  Newton-Raphson  method  can  sometimes  be  prohibitively 

delayed  if  for  some  reason  a  large  forward  bias  (V>0.8V)  is  applied  to  any 

diode  or  transistor  junction  that  has  been  represented  by  the  user  in  the 

conventional  closed  form  J  ■  I  (^^-1).  In  this  case,  the  slope  of  the 

s 

diode  curve,  given  by  m  0Ige®V,  wil1  contribute  a  very  large  term 

in  Ggg  and  consequently  in  the  Z  matrix.  The  practical  results  of  this  can 
be  more  easily  appreciated  by  consideration  of  a  one  equation  system  as 
represented  by  the  second  equation  in  Appendix  D.  A  large  derivative  caused 
by  a  highly  forward-biased  diode  leads  to  a  very  small  step  (Ax)  in  the 
Independent  variable.  Many  steps  will  therefore  be  required  to  complete 
convergence. 

If  convergence  has  occurred,  quantities  Vg,  V2,  and  Vj  are  now  known. 
There  remains  to  compute  only  capacitor  voltages  and  inductor  currents  since 
capacitor  currents  and  inductor  voltages  must  be  zero. 


2. 4. 1.2  Computing  Capacitor  Voltages 

From  the  B  matrix  in  Table  II,  the  capacitor  link  voltages  can  be 
written  in  terms  of  the  tree  branch  voltages  as 


-B14  V4  "  B15  V5  *  B17  E7 


(74) 


In  addition  the  principle  of  conservation  of  charge  permits 


-B14  C11V1+C44V4 


-B 


14 


C11  V1  (0)  +  C44  V4  (0) 


(75) 


where  (0)  and  (0)  are  initial  voltages  that  may  be  specified  by  the  user. 
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Equations  (74)  and  (75)  lead  to 


V4  *  Hc  _1  |-BU  T  CU  B15  V5  -  B14  1  C11  B17  V  BU  \l  V1  <0) 

+C44V4(0)}  (76) 

T 

where  M  -  B..  C,.  B. .  +  C.  . 

c  14  11  14  44 

Once  is  determined,  may  be  found  from  Equation  (74). 

.Note  that  Class  -4  elements  can  occur  only  if  a  capacitor  cut  set 
exists.  '  Otherwise  Equation  (76)  will  be  identically  zero  and  all  capacitor 
voltages  can  be  determined  from  Equation  (74) . 

The  use  of  Equation  (75)  permits  solution  of  a  class  of  networks  in 
which  the  final  value  of  capacitor  voltage  is  dependent  on  an  intial  value. 
Consider  the  circuit  shown  in  figure  7,  which  contains  a  capacitor  cut  set. 


Figure  7.  Capacitor  Cut  Set  Circuit 


If  both  capacitors  are  initially  uncharged,  the  final  values  of  the  capacitor 
voltages  after  the  switch  is  closed  must  be  VC1  *  1  volt,  VC2  *  9  volts.  If 
however,  capacitor  Cl  has  an  initial  charge  of  5  volts,  the  principle  of 
conservation  of  charge  (reflected  in  equation  (75))  requires  the  final  result 
to  be  VC1  ■  5.5  volts,  VC2  ■  4.5  volts. 

2.4. 1.3  Computing  Inductor  Currents 

Since  is  known, 

12  “  R22  V2  (77) 

the  inductor  tree  branch  currents  than  can  be  written  from  the  B  matrix  in 
terms  of  link  currents  as 

"*■  ^ i  Jo  ■*"  J„  +  B_,  J_  (78) 


In  addition,  the  flux  relations  around  an  inductor  loop  permit 


B36  L66 

*6  +  L33  *3  "  B36  L66  *6  (0)  +  L33  *3  (0) 

(79) 

where  I6  (0)  and  I-  (0)  are  initial  currents  that  may  be  specified  by  the 
user.  Equations  (78)  and  (79)  lead  to 

*3  ‘  V1  | 

L33  r3  (0)  +  B36  L66  X6  (0) 

“B36  L66  [ 

'»26  1  !2  +  B86  Ij8  +  B96  T  J9  +  B06  ? 

(80) 

where  ^  -  B3(. 

T 

L66  B36  +  L33 

Once  I3  is  determined,  Ig  may  be  found  from  Equation  (78).  Note  that  Class-3 
elements  can  occur  only  if  an  inductor  loop  exists.*  Otherwise,  Equation  (80) 
will  be  identically  zero  and  all  inductor  currents  can  be  determined  from 
equation  (78). 

The  use  of  Equation  (80)  permits  the  solution  of  a  class  of  networks 
in  which  the  final  values  of  inductor  currents  are  dependent  on  the  sizes  of 
the  respective  inductances.  Consider  the  circuit  shown  in  figure  8,  which 
contains  an  inductor  loop.  If  both  inductors  are  initially  relaxed,  the 
final  values  of  the  inductor  currents  after  the  switch  is  closed  must  be 
IL1  *  1  amp,  IL2  ■  9  amps. 


L2  =  1h 


Figure  8.  Inductor  Loop  Circuit 


See  paragraph  (subsequent  subsection  2. 4. 1.4)  on  network  restrictions  for 
qualification  of  this. 
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The  task  of  computing  the  capacitor  voltages  and  inductor  currents 
complete  the  initial-conditions  problem.  These  quantities  may  then  be 
transferred  to  the  transient  program  to  serve  as  initial  conditions. 

2. 4. 1.4  Restrictions 

Since  the  initial-conditions  program  is  formulated  differently  and 
for  a  different  purpose  than  is  the  transient  program,  certain  restrictions 
apply  to  the  initial-conditions  program  that  do  not  apply  to  the  transient 
program.  These  restrictions  and  the  practical  considerations  that  lead  to 
them  are  as  follows: 

1.  No  circuit  containing  a  loop  composed  entirely  of  voltage 
sources  and  inductors  will  be  accommodated.  This  situation 
would  cause  an  infinite  inductor  current  and  is  obviously  of 
no  pract  ical  importance.  The  presence  of  an  E-L  loop  is 
disclosed  by  the  condition  i  0. 

2.  No  circuit  containing  a  cut  set  composed  entirely  of  current 
sources  and  capacitors  will  be  accommodated  (see  figure  9). 
This  situation  would  invalidate  equation  (65)  and  complicate 
the  solution  process  by  requiring  that  be  carried  along  in 
the  Newton-Raphson  procedure.  This  cut  set  situation  can 
always  be  removed  by  arbitrarily  connecting  a  large  resistor 
from  node  A  to  the  ground.  Note  that  the  configuration  of 
figure  9  could  be  handled  by  the  transient  portion  of  SCEPTRE 
if  a  Newton-Raphson  solution  was  not  desired.  The  presence  of 
a  J  -  C  cut  set  is  disclosed  by  the  condition  that  either 


Figure  9.  A  Current-Source  Capacitor  Cut  Set 
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3.  The  choice  of  independent  variables  will  be  somewhat  restricted. 

No  resistor  or  inductor  current  may  be  used  as  an  independent 
variable.  Furthermore,  a  capacitor  voltage  may  be  used  as  an 
independent  variable  only  if  it  is  in  parallel  with  a  resistor 
or  current  source  (Vc  ■  Vr,  Vc  ■  Vj).  The  objective  here  is  to 
avoid  the  need  for  any  auxiliary  computation  between  passes  of 
the  iteration  procedure. 

4.  Networks  containing  only  capacitor  cut  sets  can  be  accommodated 
only  if  the  members  of  the  cut  set  are  constant. 

5.  Networks  containing  only  inductor  loops  can  be  accommodated 
only  if  the  members  of  the  loop  are  constant. 

2.4.2  MONTE  CARLO 

Monte  Carlo  calculations  in  SCEPTRE  consist  of  multiple  initial 
conditions  solutions  for  which  certain  element  values  have  been  chosen 
using  probability  distributions  specified  by  the  user.  These  variable 
elements,  representing  perhaps  production  deviations  from  nominal  values, 
cause  corresponding  variations  in  circuit  output  quantities.  The  Monte 
Carlo  function  defines  the  selection  of  variable  element  values,  controls 
the  corresponding  iterative  solutions,  monitors  output  values,  and  tabulates 
statistics  for  the  output  quantities. 

The  user  can  specify  either  the  uniform  or  Gaussian  distribution  for 
the  variable  elements.  Since  the  solution  is  calculated  for  TIME  •  0,  it 
is  independent  of  capacitive  and  inductive  elements.  Thus,  it  is  only  useful 
to  vary  resistors,  independent  current  and  voltage  sources,  and  primary  and 
secondary  current  sources  specified  by  defined  parameters. 

The  value  of  a  variable  element  is  determined  by  its  specified  mean, 
fj ,  standard  deviation,  a  ,  and  a  random  number,  picked  from  the  selected 
probability  distribution.  The  random  number  generator  provides  samples 
from  the  selected  type  of  distribution  with  zero  mean  and  unit  standard 
deviation.  Then  the  element  value  e.g.  R,  is  calculated  as 

R  «  W 

If  the  random  behavior  of  the  element  value  cannot  be  directly  defined 
in  terms  of  one  of  the  two  standard  distributions,  the  mechanism  of  defined 
parameters  may  be  used.  By  specifying  the  means  and  standard  deviations 
for  one  or  more  defined  parameters,  and  defining  the  element  value  using  an 
equation,  expression  or  table  in  terms  of  these  parameters,  the  user  can  obtain 
a  more  complicated  probabilistic  behavior  of  the  element  value. 
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As  each  new  set  of  values  Is  chosen  for  the  random  element  value, 
an  initial  conditions  calculation  is  executed.  For  each  output  quantity, 
the  sum  and  sum  of  squares  are  accumulated,  and  the  minimum  and  maximum 
observed  over  all  Iterations  are  tabulated.  The  mean  and  standard  deviation 
for  each  output  X  are  then  computed  as 


and 


M  - 

x 


£  X 

i  i 


(81) 


1/2 


(82) 


where  is  the  mean,  is  the  standard  deviation,  and  the  sums  are  taken 
from  1  to  N,  where  N  is  the  number  of  Monte  Carlo  iterations. 

2.4.3  SENSITIVITY 

The  Sensitivity  option  gives  the  user  the  partial  derivative  of  a 
network  function  with  respect  to  a  list  of  network  independent  variables. 
This  option  gives  the  user  the  "un- normalized  sensitivity",  i.e. 


(83) 


The  normalized  sensitivity,  or  the  percentage  change  in  a  function  with 
respect  to  the  percentage  change  in  an  independent  variable,  may  be  obtained, 
if  desired,  by  using  a  defined  parameter.  The  normalized  sensitivity  is 


To  calculate  sensitivity,  an  initial  conditions  run  is  made  on  the 
original  network,  and  appropriate  network  element  currents  and  voltages  are 
saved.  An  adjoint  network  as  described  in  Appendix  F  is  generated.  (Ref.  10, 
11.)  For  each  network  function  whose  sensitivity  is  requested,  an  initial  con¬ 
ditions  run  is  made  on  the  adjoint  network  with  a  proper  source,  and  the  element 
currents  and  voltages  in  the  adjoint  run  are  saved.  The  proper  product  of  the 
original  network  values  and  the  adjoint  network  values  (see  Appendix  F)  gives 
the  required  sensitivity.  Table  III  gives  the  sensitivity  (superscript  "a" 
indicates  adjoint  values). 
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TABLE  III 


SENSITIVITY  COMPUTATIONS 


INDEPENDENT  VARIABLE 

SENSITIVITY 

Resistor  (R) 

(Current  through  the  resistor  in  original 

IC  run)  x  (Current  through  the  same 
resistor  in  adjoint  IC  run),  i.e., 

<IR)  (IRa) 

Independent  Voltage 

Source  (E^) 

Negative  current  through  the  voltage 
source  in  adjoint  IC  run  i.e.,  IE| 

Independent  Current 

Source  (Jg) 

Negative  voltage  across  the  current  source 
in  adjoint  IC  run  i.e.,  -VJj* 

O 

Factor  in  Secondary 

Current  Source  given  by  Defined 
Parameter  (P)  where  JQ  ■  P*J^ 

(The  primary  current  in  original  IC  run) 
x  (Negative  voltage  across  the  secondary 
source  in  adjoint  IC  run)  i.e.,  (JQ) 

(~VJg)  9 

_ _ _ _ _ _ _ _  ,  .  I.. 
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2.4.4  WORST  CASE 


A  Worst-Case  calculation  extrapolates  linearly  along  the  gradient  vector 
of  an  objective  function,  F,  to  reach  a  function  minimum  and  maximum  at  one  or 
more  independent  variable  bounds.  The  gradient  vector  VF  is  computed*  for  nominal 
circuit  conditions,  and  it  is  assumed  that  F  depends  linearly  on  the  N  independent 
variables.  The  gradient  indicates  the  direction  in  which  the  function  grows  most 
rapidly,  and  the  assumed  linearity  allows  extrapolation  along  this  vector  to  reach 
the  minimum  and  maximum. 

The  maximum  of  F  is  then  obtained  by  extrapolating  the  independent  variable 
values  X  as  far  as  possible  along  the  gradient  vector: 


^  «  Y  +  XVF, 


\>0 


(85) 


Here,  Xjj  is  the  independent  variable  vector  producing  the  high  value  of  F,  Y 
is  the  vector  of  independent  variables  at  the  nominal  point,  and  \  is  taken  as 
large  as  possible  with  Xh  remaining  within  bounds. 

Under  the  same  assumption,  the  minimum  of  F  is  obtained  by  extrapolation 
along  the  negative  of  the  gradient  vector: 


XL  =  Y  -A/VF,  /J>0 


(86) 


Xl  is  the  vector  of  independent  variables  at  the  point  producing  the  low 
value  of  F,  and  u  is  chosen  as  large  as  possible,  keeping  the  independent 
variables  within  bounds. 

The  independent  variables  are  specified  with  lower  and  upper  bounds 


P  £  X  £  i  -  1  ...N 


See  Appendix?  for  a  discussion  of  the  Adjoint  method,  used  for  computing 
the  gradient  vector. 
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If  the  i —  component  of  the  gradient  vector 
increases  as  increases  from  Y^  toward  Q^. 

F  increases  as  decreases  from  Yi  toward  P^. 


—rr  is  positive,  then  F 

<*xi  9F 

If*  however,  -gj 


is  negative, 


as 


The  general  equation  of  a  line  in  N  dimensional  space  may  be  written 


xrYi 


x2-y2 


Vyn 


M. 


where  (X^,  i«l...  N)  and  (Y^,  i*l...N)  are  two  points  on  the  line  and  Jj  is 
a  unit  vector  in  the  direction  of  the  line. 


The  equation  can  be  restated 


Vli 


■  constant 


rth 


(87) 


When  the  J—  component  takes  on  its  lower  bounding  value,  we  have 


VYi 


P  -Y 

llJj- 


,  i  -  1...N 


or. 


Xi  *  ^  <W  +  Yi 


The  length  of  the  line  segment  is  given  by 

.2 


s  -  I  <w2  - 


v 

^  1 


(W 


since  /j  is  a  unit  vector. 
Thus , 


j 


P  -Y 

1-1 
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the  distance  along  the  gradient  to  the  lower  bounding  surface  of  the  J — 
independent  variable. 


Similarly, 


2£i 


» 


the  distance  to  the  corresponding  upper  boundary  surface. 


For  the  high  limit  calculation,  a  search  is  made  through  all  Dp^  for 

dF  o  p 

which  *s  positive,  and  all  Dq^  for  which  *s  ne8ative,  until  the 

smallest  number  from  this  set  is  found.  (Note  that  F  is  independent  of  any 

9f 


X^  for  which 


“5X1 


■  0) .  This  minimum  value ,  D. 


HIGH’ 


represents  the  largest 


distance  X  can  move  in  the  direction  of  the  gradient  vector  before  it  reaches 
a  boundary.  Xjjjqjj  is  the  independent  variable  which  established  the  allowable 
displacement  of  X  associated  with  DUT(_,. 

HIGH 


A  similar  search  through  all  Dp^  for  which 


ax 


is  negative,  and  all 


Dq^  for  which 


which  is  the 


is  positive,  produces  another  minimum,  DL0W* 

largest  possible  displacement  of  X  from  Y  along  the  negative  of  the  gradient 
vector  before  striking  a  boundary.  X^qw  is  the  independent  variable  which 
established  the  allowable  displacement  associated  with  DL0W- 


XHIGH  a°d  xlow  may  or  may  not  be  the  same  variable.  Figure  10 
illustrates  a  two-dimensional  case  where  the  high  limit  is  determined  by 
the  upper  bound  of  one  variable  (X2) ,  and  the  low  limit  is  determined  by 
the  upper  bound  of  another  (Xjj . 


When  the  distance  D„T  has  been  determined,  Equation  (87)  can  be  used 
to  define  all  independent  variable  values,  Xjjjqjj,  at  the  high  limit. 


or 


VYi 


HIGH 


*H 


Y  +  D 


HIGH 


(88) 
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\ 


\ 


HIGH  LIMIT 
DETERMINED  BY  Q, 

XHIGH  IS  X2 


LOW  LIMIT 


k 


\ 


Q, 


Figure  10.  Two-Dimensional  Example  of  Worst-Case  Limits 


Note  that  this  is  exactly  in  the  form  of  Equation  85,  since  y  is  pro¬ 
portional  to  VF  and  D^j^  >0. 

In  a  similar  way,  we  find  the  independent  variables  at  the  low  limit 

to  be 


V7-dlow  "  (85> 

as  required  by  Equation  (86). 

When  Xy  and  have  been  found,  computations  of  objective  function 
value  and  gradient  at  each  of  the  two  points  are  carried  out.  Inspection 
of  the  gradient  vector  at  these  points  can  give  information  on  whether  the 
true  extreme  values  occur  within  the  bounded  area  or  if  they  are  external. 

2.4.5  OPTIMIZATION 

The  Optimization  option  finds  the  local  minimum  of  an  objective  function 
with  respect  to  a  list  of  bounded  independent  variables.  The  local  minimum  is 
the  point  in  the  independent  variable  space  at  which  each  component  of  the 
gradient  either  is  zero,  or  is  negative  (function  decreasing)  with  its 
corresponding  independent  variable  at  a  bound. 

The  Optimization,  or  more  accurately  minimization,  function  is  based 
on  the  Davidon  algorithm  (Ref.  13).  As  in  most  such  processes^  the  gradient  is 
computed  for  some  set  of  values  of  the  Independent  variables,  Y,  within  all 
boundaries,  and  the  independent  variable  vector  is  displaced  along  the  gradient. 
At  this  new  point,  the  gradient  vector  is  again  computed,  and  the  cycle  repeats 
until  the  error  (discussed  below)  is  less  than  some  criterion  value. 

The  adjoint  method  of  computing  tne  gradient  obviates  the  use  of  a 
direct  search  method  with  its  inherent  numerical  problems.  Furthermore, 
since  the  adjoint  procedure  determines  the  gradient  with  only  two  IC 
solutions  regardless  of  the  number  of  independent  variables,  it  has  a 
speed  advantage  over  direct  search  methods  and  other  gradient  methods. 

This  speed  advantage  is  enhanced  because  the  auxiliary  IC  solution  (for 
the  adjoint  network)  is  strictly  linear  and  will  generally  converge  quickly. 
Thus  the  combination  of  rapid  gradient  calculation  and  efficient  search  for 
the  function  minimum  provide  the  SCEPTRE  user  with  an  extremely  versatile 
tool. 


The  Davidon  method  is  a  synthesis  of  two  techniques:  The  Steepest 
Descent  Method  and  the  Newton-Raphson  method  applied  in  the  Initial  Conditions 
calculation  (subsection  2.4.1).  The  Newton-Raphson  method  is  based  on  the 
truncated  Taylor  expansion 

F  (X)  =  F  (Y)  +  h  '  V  F  (Y)  +  1/2  h  ‘G  (Y)  h 

where  h  is  the  displacement  vector  (X  -  Y)  and  G  is  the  Hessian  matrix  of 
second  partial  derivatives  of  the  objective  function. 

d  2f  . 

Gij  6Xi  aXj 
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Differentiation  of  this  approximation  to  F(X)  leads  to  the  expression  for  the 
gradient . 

VF(X)  -  VFCY)  +  G(Y)  "h. 

At  the  minimum  value  X,  TF(X)  ■  0,  so  that  the  displacement  h  is  given  by 
TT  -  -G"\7)  VF(Y). 

In  the  vicinity  of  the  minimum,  it  is  assumed  that  F  is  essentially  a  quadratic 
function  of  X,  so  that  the  point  at  which  G  is  evaluated  is  not  critical. 

The  Steepest  Descent  step  is  determined  completely  by  the  gradient 
vector,  and  thus  has  the  advantage  of  simplicity.  However,  the  process  is 
very  inefficient,  even  for  functions  with  quadratic  dependence  on  the 
independent  variables,  and  is  not  generally  recommended  except  for  initiating 
a  minimization  calculation. 


Because  the  Davidon  technique  smoothly  changes  from  the  Steepest  Descent 
to  the  Newton-Raphson  iteration,  maintaining  suitable  conditions  on  h  at  all 
times,  it  is  a  very  effective  procedure  (Ref.  14). 

The  Davidon  process  (Ref.  15)  starts  with  an  initial  approximation  H0  to  G 
which  is  generally  the  unit  matrix  and  is  thus  positive  definite  (since  the 
eigenvalues  all  have  the  value  unity).  The  iterations  produce  an  increasingly 
better  approximation  H  to  G“1  without  matrix  inversion,  while  continuing  to 
enforce  the  condition  of  positive  definiteness  which  will  be  required  at  the 
minimum. 


The  iteration  for  the  Davidon  method  is  given  by 
1(1+1)  =  X(i)  +  h^  Tf^ 

where 

X^  is  the  vector  of  independent  variables; 

h^  is  the  step  size,  chosen  by  an  auxiliary  linear  search  process; 
is  the  approximation  to  G  ^ 

and  the  superscript  identifies  the  iteration  at  which  the  function  is  evaluated. 
The  matrix  H  is  updated  as  follows.  Define  the  vectors 

Z  — h(1)  H(i)  VF(i) 


ZT=  vf^1+1^  -  VF^ 

nd+i)  _  H(i)  +  A(i)  +  g(i) 
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where  Che  elements  of  matrices  A 


(i) 


and  B 


U> 


are  respectively  given  by 


■ 

Z=A 

1  •  7i 

‘S?  ■ 

(H(1Wl(H(1V)k 

7T*h  77 

The  matrix  A  provides  convergence 
positive  definiteness  of  H  at  each  step. 

,-l 


while  B  maintains  the 


The  implementation  in  SCEPTRE  is  based  on  the  original  formulation 
produced  by  the  Argonne  National  Laboratory  (Ref.  16).  However,  the  original 
code  was  inappropriate  for  constrained  optimization  and  so  the  transformation 
of  Independent  variables  suggested  by  Box  (Ref.  14)  was  incorporated  to  pro¬ 
perly  take  care  of  the  Inequality  constraints  of  the  SCEPTRE  independent 
variables.  This  transformation  may  be  stated  in  the  following  way:  Suppose 
the  ith  independent  variable  is  restricted  so  that 

pi*xi*Qi 

Then  the  transformation  of  variables 


X,-  (2£i)_  (!rt) 


cos  Y, 


is  introduced.  The  corresponding  derivative  is  given  by 


Note  that 


sin 


0 


TT 


and 


at  both  extremes 


The  gradient  components  with  respect  to  the  new  variables  are  given  by 

9F  ^i  _£F 
"STi  "  dY±  ax1 
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Thus,  the  effect  of  the  transformation  Is  to  introduce  zeros  of  the  gradient 
along  the  boundary,  making  these  points  acceptable  to  the  algorithm.  No  such 
zeros  are  introduced  within  the  boundaries,  and  the  algorithm  will  preferentially 
converge  to  an  internal  minimum  if  one  exists. 

Another  feature  is  the  removal  of  a  restriction  present  in  the  original 
formulation,  which  required  the  minimum  value  of  the  objective  function  to  be 
positive.  A  floating  reference  value  has  been  introduced  which  is  smaller 
than  the  estimate  of  the  minimum  at  any  stage  of  the  iteration  and  which  is 
automatically  updated  as  required.  The  user  may  specify  a  reference  value 
known  to  be  smaller  than  the  actual  minimum.  If  he  does  so,  the  process 
will  probably  converge  more  quickly,  although  a  poor  choice  may  increase  the 
number  of  Iterations  required. 


Figure  11.  Change  in  AF  Due  to  Change  in  Independent  Variable 

The  Iteration  stops  when  the  change  in  function  value  expected  from 
the  next  increment  to  the  Independent  variable  vector  is  less  than  the 
specified  tolerance. 

Thus ,  in  figure  LL,  if  the  anticipated  change  A.F  due  to  changing  the 
independent  variable  to  X  (1+^) ,  computed  at  point  X(i^  using  the  supposed 
quadratic  dependence  of  F  cn  X,  is  less  than  the  tolerance,  the  value  X^+^^ 
will  be  reported  out  as  the  location  of  the  minimum,  but  no  further  isolation 
will  be  done. 
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2.5  AC  SOLUTIONS 

This  section  describes  the  AC  analysis  by  SCEPTRE  for  the  following 
situations: 

Case  (1)  Circuits  Containing  Independent  Source  Only 
X  -  AX  +  GV,  where  VT  -  [E?  Jg] 

Case  (2)  Circuits  Containing  Linearly  Dependent  Sources 

’  T 

X  -  AX  +  GV  +  Hu,  where  u  -  [EY  JX] 

Case  (3)  Circuits  Requiring  Time  Derivatives  of  Independent  Sources 

X  =  AX  +  GV  +  QV 


Case  (4)  Circuits  Containing  Linearly  Dependent  Sources  and  Requiring 
Time  Derivatives  of  Independent  Sources 


X  *  AX  +  GV  +  QV  +  Hu 


Circuits  containing  linearly  dependent  sources  and  requiring  time  derivatives 
of  both  independent  and  dependent  sources  are  excluded  from  both  transient 
and  AC  analysis. 

2.5.1  INDEPENDENT  SOURCES  ONLY 


2. 5. 1.1  State  Variable  Formulation 


Consider  the  frequency-domain  analysis  of  linear  time-invariant 
systems  governed  by  the  equations: 

X  =  AX  +  GV  (90) 

Y  «  CX  +  DV  (91) 

where 

A  is  an  (N  x  N)  matrix  of  real  coefficients  describing  the  network. 

G  is  an  (N  x  NS)  matrix  of  real  coefficients  relating  the 
externally  applied  sources  to  the  network. 

N  is  the  number  of  state  variables. 

NS  is  the  number  of  externally  applied  sources. 

NE  is  the  total  number  of  elements  in  the  network. 

NR  is  the  number  of  outputs  requested,  and  NR<2  *  NE. 
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C  is  an  (NR  x  N)  matrix  of  real  coefficients  relating  the 
state  variables  to  the  outputs. 

D  is  an  (NR  x  NS)  matrix  of  real  coefficients  relating  the  inputs 
to  the  outputs. 

X  is  an  (N  x  1)  vector  composed  of  state  variables. 

V  is  an  (NS  x  1)  vector  composed  of  externally  applied  sources. 

V  is  an  (NR  x  1)  vector  composed  of  requested  outputs. 

A  method  which  will  permit  multiple  inputs,  i.e.,  multiple  penetra¬ 
tions,  with  frequency  dependent  magnitude  and  phase  variation.  For  steady- 
state  AC  analysis,  this  requirement  on  each  individual  driving  function, 
lSk^NS,  may  be  indicated  as  follows! 

Vk(t,u)  -  ^  (u :)  e  at  (92) 

where 

a'  is  frequency  in  radians. 

Q  is  a  complex  quantity  given  by  a  ■  ju,’. 

K.  (u.)  is,  in  general,  a  complex  function  of  frequency. 

Writing  Equation  (92)  in  vector  notation  we  have 

V  -  Ke  Qt  (93) 

where  the  (NS  x  1)  column  vector  V  is  written  as  NS  terms  of  the  complex 
frequency  variable  Q  with  each  term  having  an  independent  magnitude  and 
phase  variation  given  by  K. 

The  solution  of  Equation  (90)  is  of  the  form 

X  «  X  e  Qt  (94) 

o 

Using  Equations  (93)  and  (94)  we  can  rewrite  equation  (90) 
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and 


or 


OX  -  AX  +  GK 
o  o 


(a I  -  A)Xq  -  GK 
Solving  for  Xq  we  get 

Xq  -  (01  -  A)-1GK  (95) 

The  system  matrix  A  can  be  reduced  by  similarity  transformations  to  a 
diagonal  matrix  A.  This  is  done  by  the  EIGENP  routine(Ref.  17,  18).  The 
process  Is 

A  -  SAS-1 

where 

\  Is  an  (N  x  N)  diagonal  matrix  composed  of  the  complex  system 
eigenvalues,  and 

S  is  a  complex  (N  x  N)  matrix  whose  columns  are  the  system 

eigenvectors.  This  matrix  is  referred  to  as  the  modal  matrix. 

The  substitution  of  (96)  into  (95)  leads  to 

X  -  (SQS_1  -  S  \S_1)-1GK 
o 

X  -  S(Q  I  -  A)"1s'1gk 

O 

or 

Xo  -  S(j  u)I  -  \)-1S-1GK(uO  (97) 

Equation  (97)  represents  the  AC  solution  to  the  state  variable  equation 
given  in  (90) . 

Equation  (91)  represents  the  non-state  variable  quantities.  The 
AC  solutions  for  these  voltages  and  currents  are  consistent  with  the 
transient  equations  given  in  subsection  2.3.1,  with  the  understanding  that 
their  implementation  requires  complex  arithmetic.  The  simpler  equations 
below  are  substituted  for  four  quantities. 

1*1  <JL>C  V 
L1  J  U11  1 

h  ■ 

V3  *  (L36!6  +  L33V 

v6  ■  <L36:3  +  L66V 
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2. 5. 1.2  A  and  G  Matrix  Calculations 


SCEPTRE  defines  the  state  variables  as  the  link  inductor  currents, 
I3,  and  the  capacitor  branch  voltages,  V4.  SCEPTRE  also  permits  both 
voltage  sources,  E7,  and  current  sources,  Jg.  This  leads  to  the  following 
definition  for  X  and  V. 


Rewriting  and  partioning  Equation  (90)  in  an  alternate  matrix  form 

produces : 


as  separate  equations  we  have 


X3  “  AH  *3  +  A12  V4  +  G11  E7  +  G12  J8 


V4  "  A21  *3  +  A22  V4  +  G21  E7  +  °22  J8 

The  transient  equations  of  subsection  2.3.1  (Ref.  17,  18)  were  used 
in  the  derivation  of  the  necessary  submatrices. 

The  submatrices  A^.,  A A^  and  A^  ^or  t*ie  case  Involving  lumped, 
passive,  time  invariant  ?,  L  and  C  elements  are  given  by  Equations  (98), 

(99),  (100),  and  (101). 

An  -  (-MLI)*B35*R55)  *  |"b35  -  B^  *  (MRI)  *  B25*R55*B35 1  (98) 


35*R55>  *  [B35  -  B25  *  (MRI>  *  B25*R55*BI5] 
A12  -  (-MU)  *[b34  -  B35*R55*B^MMRI)*B24J 

[b34  *  (msi)*b14*s44  -  il  *  ^4*(KRI)*B25*I 


A21  ■  <C44»*  Bl«  *  <HSI)*BU*SU 


B24*(MRI)*B25*R55*B35  -  B 

(100) 


A22  -  (c44x)  *  B14*(MSI)*B14*S44 


-  1]  *[b24  *b241  (101) 
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where 
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MLI  -  I  L33  +  Lj6*bJ6  +  B36*L6j  +  B36.L66*B^6 


[■ 


]" 


MR  I  - 


MSI  - 


R  +  B  *R  *B 
r22  d25  55  b25 


S11  +  WB14 


-1 


-1 


The  submatrices  Gj^,  Gj^,  and  ^21»  G22»  ^or  t^ie  case  involving 
independent  voltage  and/or  current  sources,  are  given  by  Equations  (102) 
through  (105). 
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(-MLI)  *  |^B37  -  B35*R55*B25*(MRI)*B27J  (102) 

(-MLI)  *  |b35*R55*(bJ5  -  B25*(MRI)*B25*R55*Bg5J  (103) 

(C44r)*  [(B14*(MSI)*B14*S44  "  1  > *B24* <MRI) *B27 ]  U04) 

(C44X)  *  [(BU*(MSI)*B14*S44  -  I)*(B^(MRI)*B25*R55*B^5  -  B^J 

(105) 


2.5.2  CIRCUIVS  CONTAINING  LINEARLY  DEPENDENT  SOURCES 
2. 5. 2.1  Reduction  to  Canonical  Form 


SCEPTRE  permits  the  user  to  define  voltage  and  current  sources 
that  are  linear  functions  of  resistor  voltages  and  currents.  These  are 
type  EY  and  JX  sources.  Let  u  be  a  vector  composed  or.  these  sources 
T 

given  by,  u  =  [EY  JX] .  Then  a  circuit  containing  thise  sources  can  be 
characterized  by 

X  =  AX  +  GV  +  Hu  (106) 

where  H  is  an  (N  X  ND)  matrix  of  real  coefficients  relating  the  dependent 
sources  to  the  state  variables. 

ND  is  the  number  of  linearly  dependent  sources.  The  problem  is  to 
reduce  Equation  (106)  to  canonical  form,  namely 

X  -  (AEXT)X  +  (GEXT)V 
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where  AEXT  is  an  (NXN)  matrix  of  real  coefficients  describing  the  extended 
network.  GEXT  is  an  (N  X  NS)  matrix  of  real  coefficients  relating  the 
external  sources  to  the  extended  network.  Equation  (106)  may  be  rewritten 
in  the  following  way: 
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By  definition. 


(107) 


EY 


*  klV2 


k2V5 


and 


(108) 


JX  -  k3I2  +  k4I5  (109) 

where  k^,  k2»  k^  and  k^  are  matrices  relating  each  dependent  source  to 
the  appropriate  resistor  voltage  or  current  on  which  it  depends. 


Since  V.  =  R0„I.  and  I,  =  G,,V,  substitution  into  (108)  and  (109) 
yields:  2  22  2  5  55  5 


EY 

JX 


klR22 


(110) 


Substituting  (110)  into  (1fl7)  gives 
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(111) 

T 

The  vector  [Y]  =  [1^^]  can  be  expressed  in  terms  of  state  variables 

and  independent  sources  by  observing  that  Y^  is  a  subset  of  the  non-state 
variable  outputs  Y,  where  Y  =  CX  +  DV.  That  is 


c'x  +  d'v 


(112) 
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However , 


subsection  2. 3. A  gives 
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(113) 


(114; 


Using  (114),  (113)  and  (112)  we  can  write  Equation  (111)  in  terms 
of  state  variables  and  independent  sources 
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Collecting  terms  we  have 


where  Equation  (115  represents  a  reduction  of  Equation  (106)  to  canonical 
form  as  required  with 


2.5.2. 2  H  Matrix  Calculation 

Eauation  (115)  is  complete  except  for  a  description  of  the  H  matrix 
which  is  necessary  for  computational  purposes. 

It  is  easily  shown  that  the  H  matrix  is  quite  similar  in  structure 
to  the  G  matrix  by  observing  the  equivalence  in  the  following  equation 


This  indicates  that  the  G  and  H  submatrices  differ  only  by  post¬ 
multiplications  by  terms  which  relate  the  circuit  connectivity  to  the 
appropriate  sources. 
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Thus 
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2.5.3  CIRCUITS  REQUIRING  TIME  DERIVATIVES  OF  INDEPENDENT  SOURCES 
2. 5. 3.1  Derivation  of  the  AC  Solution 


If  the  circuit  topology  is  such  that  a  source  time  derivative  is 
required  in  order  to  perform  a  transient  analvsis*,  then  this  same  require¬ 
ment  is  also  imposed  on  an  AC  analysis.  However,  there  is  one  important 
difference.  In  the  AC  analysis  program,  unlike  the  transient  program,  the 
user  is  not  required  to  supply  the  derivative.  The  AC  program  automatically 
supplies  the  needed  information  and  proceeds. 

The  derivatives  are  treated  internally  as  follows: 


*  Independent  voltage  sources 
Independent  Current  Sources 


together  with  a  capacitor  tie-set,  or 
together  with  an  inductor  cut-set. 
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Let  this  situation  be  characterized  by 


•  • 

X  -  AX  +  GV  +  QV  where, 

Vk(t,u')  -  Kk(^)eJa’t 

Then  the  solution  of  Equation  (116)  is  of  the  form 
1  u.'  t 

X  -  X  eJ 
o 


(116) 


Thus, 

-j-  (X  eJa,t)  *  A(Xeju't)  +  G(KeJu't)  +  Q~(KeJut) 

dt  O  O  at 

Ja’X  -  AX  eJw,  t  +  GKeJ  C  +  j  u? QKeJ  a'  l 

o  o 


and 

(ja'I-A)XQ  =  GK  +  la'QK 
Xq  =  (ja'I  -  A)”1  (GK  +  j^QK) 
Substituting  our  similarity  transformation  gives 


Xq  =  S  ( j  u,'  I  -  N)'1  S-1(GK  +  ja'QK) 

X  =  S(jo,'I  -  \)_1  [(S_1GK)  +  j  u,'  (S_1QK)  ]  (117) 

o 

Equation  (117)  represents  the  AC  solution  to  the  state  variable  Equation 
given  in  (116),  The  second  term  in  the  brackets  represents  the  additional 
computation  required,  at  each  frequency,  due  to  the  presence  of  all  source 
derivatives. 


2. 5. 3. 2  Q  Matrix  Calculation 


This  matrix  was  also  derived,  using  the  differential  equations  of 
Section  2.3.1.  It  is  listed  here  for  continuity. 
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where 


2.5.4  CIRCUITS  CONTAINING  LINEARLY  DEPENDENT  SOURCES  AND  REQUIRING  TIME 
DERIVATIVES  OF  INDEPENDENT  SOURCES 

This  case  presents  no  difficulty.  It  is  simply  a  combination  of 
Cases  (2)  and  (3). 
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SECTION  III 


SYSTEM  OPERATION 


3.1  INTRODUCTION 


The  SCEPTRE  System  is  a  FORTRAN  IV  computer  program  written  for  the  IBM 
7090/94  and  360  Data  Processing  Systems.  It  consists  of  two  major  phases. 

The  first,  called  the  Program  Generator,  creates  (on  a  disk  or  tape  data 
file*)  another  FORTRAN  IV  computer  program  containing  circuit  equations  for 
electrical  networks.  SCEPTRE  does  this  automatically,  using  the  input  data 
describing  the  circuit  to  be  analyzed.  The  second  phase,  the  Circuit  Solution 
Executive,  computes  the  circuit  response  by  solving  these  equations  generated 
in  the  FORTRAN  IV  program. 

Operation  of  the  SCEPTRE  System  depends  upon  the  Monitor  System,  IBSYS 
(operating  system  OS/360)  which  controls  its  execution.  Under  control  of 
IBSYS  (OS/360),  SCEPTRE  executes  in  two  phases  as  separate  job  steps  that 
are  loaded  and  executed  sequentially.  Prior  to  loading,  however,  IBSYS  (OS/ 
360)  performs  any  required  FORTRAN  compilation.  Programs  to  be  compiled  do 
not  have  to  reside  on  the  standard  input  file  upon  which  the  input  data  are 
stored.  Instead,  the  system  can  be  instructed,  via  the  job  control  language, 
JCL  to  compile  the  load  programs  from  an  alternate  input  file.  SCEPTRE  uses 
this  feature  of  IBSYS  (OS/360)  as  a  means  of  linking  its  two  job  steps.  This 
is  illustrated  in  the  System  Flow  Diagram,  Figure  12. 

3.2  PROGRAM  GENERATOR 


The  Program  Generator  is  an  executive  program  that  controls  the  inputting 
of  circuit  description  data,  generation  of  a  FORTRAN  IV  subprogram  for  calc¬ 
ulating  circuit  response,  generation  of  circuit  parameter  data,  storage  of 
circuit  models  on  library  files,  restart  of  discontinued  runs,  and  re-outputting 
of  computed  results.  Each  of  these  six  program  tasks  as  well  as  the  Program 
Generator  (EXEC1)  are  described  in  the  flow  diagrams,  figures  13  through  18. 

3.2.1  CIRCUIT  DESCRIPTION  PROCESSOR 

The  SCEPTRE  circuit  description  language  is  used  to  describe  electrical 
networks.  This  application-oriented  language  is  powerful,  easy  to  learn  and 
use,  and  nearly  format  free.  With  it,  circuits  composed  of  fundamental 
circuit  elements  and  prestored  circuit  models  can  be  described.  The  types 
of  fundamental  electrical  characteristics  allowed  are  resistance,  capacitance, 
inductance,  current  and  voltage  sources,  and  mutual  inductance.  Using  the 
same  components  and  circuit  description  language,  equivalent  circuit  models 
of  devices  such  as  transistors  and  diodes  can  be  described  and  stored  on  a 
library  file  for  future  use.  When  a  stored  model  is  referenced  in  a  circuit 
description,  it  is  located  on  the  specified  library  file  and  its  elements 
appropriately  substituted  into  the  circuit  being  described. 

*This  term  "file"  in  this  section  refers  to  disk  or  tape  file  for  the 
S/360  and  to  a  tape  file  for  the  7090/94. 
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Figure  15.  Model  Editor  (MDEDIT)  Flow  Diagram 
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Figure  18.  Circuit  Solution  Executive  (EXEC  2)  Flow  Diagram 


Often,  when  a  model  is  called  out  in  describing  a  circuit,  the  desired 
parameter  values  are  different  from  those  originally  specified.  Rather  than 
permanently  change  the  stored  model,  the  SCEPTRE  circuit  description  language 
allows  model  parameter  values  to  be  changed  when  model  element  substitution 
takes  place. 

3.2.2  MODEL  EDITOR 

Using  the  SCEPTRE  circuit  description  language,  circuit  models  const- 
ucted  by  the  user  can  be  stored  on  a  model  library  file  by  the  Model  Editor 
program.  One  permanent  library  file  and  one  temporary  library  file  are 
provided  for  storing  models.  Circuit  models  may  be  an  arbitrary  n-terminal 
configuration  of  the  allowed  fundamental  circuit  elements.  Any  model  so 
described  can  be  stored  on  either  library  file  using  the  Model  Editor  Program. 

3.2.3  CIRCUIT  EQUATION  GENERATOR 

After  the  description  of  the  circuit  to  be  analyzed  has  been  reconstructed 
in  memory  of  the  Circuit  Description  Processor,  the  FORTRAN  IV  sub-program  called 
SIMUL8  is  created.  It  contains  D-C  steady-state,  transient,  and/or  AC  solution 
equations,  depending  on  the  type  of  analysis  requested  by  the  user.  These 
equations  are  based  on  the  formulation  presented  in  Section  II  of  this  report. 

The  SIMUL8  program  is  written  on  an  output  file  (PROGRAM  SAVE  TAPE)  and  stored 
until  the  second  phase  of  SCEPTRE  operation,  when  it  will  be  compiled  and 
executed . 

3.2.4  DATA  GENERATOR  AND  RERUN  PROCESSOR 

Essentially  the  SIMUL8  program  contains  only  the  circuit  equations  for 
the  network  under  investigation.  The  parameter  or  component  values  of  the  net¬ 
work  are  stored  separately  as  input  data  to  the  SIMUL8  program.  The  Data  Gen¬ 
erator  program  organizes  and  stores  the  circuit  data  on  the  PROGRAM  SAVE  TAPE. 

Once  a  circuit  has  been  described  to  SCEPTRE,  multiple  or  repeated 
circuit  solutions  can  be  run  by  changing  parameter  values  between  runs.  The 
circuit  description  language  is  used  to  specify  the  number  of  repeated  runs, 
each  known  as  a  rerun,  and  the  changes  in  parameter  values  desired  for  each 
rerun.  The  Rerun  Processor  interprets  and  processes  this  information.  The 
Data  Generator  then  creates  and  stores  one  block  of  circuit  data  for  each 
rerun  on  the  PROGRAM  SAVE  TAPE.  Since  only  parameter  values  are  changed  between 
reruns  (i.e.,  no  topological  changes),  the  SIMUL8  program  containing  the  original 
circuit  equations  is  simply  re-executed  during  the  solution  phase  using  the  data 
created  for  each  rerun. 
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3.2.5  CONTINUE  PROCESSOR 


The  Continue  Processor  allows  previously  terminated  transient  solution 
runs  to  be  restarted  and  continued.  In  addition,  it  is  necessary  that  the 
PROGRAM  SAVE  TAPE  from  the  original  transient  solution  run  be  available. 

This  file  contains,  in  addition  to  the  SIMUL8  program,  all  of  the  parameter 
values  and  data  required  to  continue  the  transient  solution.  The  Continue 
Processor  also  allows  certain  run  control  parameters  to  be  changed  and  used 
for  the  continuation  of  the  original  transient  run. 

3.2.6  RE-OUTPUT  PROCESSOR 

The  SCEPTRE  user  may  wish  to  re-create  both  the  printed  and  plotted 
outputs  for  a  particular  transient  analysis  run.  This  can  be  done,  providing 
the  OUTPUT  SAVE  TAPE  is  saved  from  the  original  solution  run,  by  using  the 
Re-Output  Processor.  Several  changes  in  the  output  can  be  made  during  a  Re- 
Output  run.  For  example: 

-  Output  labels  can  be  changed 

-  New  quantities  can  be  plotted 

_  The  order  in  which  output  quantities  are  printed  out  and  plotted 
can  be  changed 

-  Printing  and  plotting  of  output  quantities  can  be  suppressed. 

3.3  SOLUTION  EXECUTIVE 


In  the  second  phase  of  SCEPTRE  operation,  the  FORTRAN  IV  sub-program 
SIMUL8  is  compiled  and  executed.  During  this  job  step,  IBSYS  (OS/360),  is 
instructed  to  read  the  SIMUL8  program  from  the  PROGRAM  SAVE  TAPE,  compile  it, 
and  then  link  edit  it  along  with  the  other  required  programs. 

Execution  commences,  as  shown  in  figure  18  with  the  reading  of  the  run 
control  parameters  from  the  PROGRAM  SAVE  TAPE.  Then  the  SIMUL8  program  is 
called.  It  reads  the  remaining  circuit  parameter  values  stored  on  the  PROGRAM 
SAVE  TAPE  and  then  enters  the  circuit  solution  equations.  The  circuit  equations 
are  solved,  thus  computing  the  state-variable  derivatives.  Calling  the  selected 
numerical  integration  routine  then  produces  new  values  of  the  state  variables  ior 
the  next  time  step.  At  the  conclusion  of  each  successful  integration  step,  the 
requested  output  quantities  are  buffered  in  memory.  When  the  buffer  is  full, 
it  is  written  on  the  OUTPUT  SAVE  TAPE  as  binary  data.  After  the  circuit  solution 
is  complete,  control  is  returned  to  the  Solution  Executive  Program,  whereupon  the 
contents  of  the  OUTPUT  SAVE  TAPE  (computed  results)  are  re-formatted  in  lists  and 
graphs  and  stored  on  the  SYSTEM  OUTPUT  TAPE  for  peripheral  processing. 

If  any  reruns  were  requested,  control  is  then  returned  to  SIMUL8  for  re- 
execution  of  the  solution  p'*ase  and  subsequent  outputting.  This  recycling  is 
continued  until  all  circuit  reruns  have  been  processed. 
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At  the  conclusion  of  each  transient  solution,  the  critical  parameter 
values  and  data  are  stored  on  the  PROGRAM  SAVE  TAPE.  Thus,  if  the  file  is 
saved  at  the  end  of  the  run,  the  solution  can  be  continued  at  some  future  time. 
Use  of  the  Continue  Processor  allows  previously  terminated  solution  runs  to  be 
restarted  and  continued  with  changes  in  the  run  control  parameters. 

Transient  solution  runs  may  be  terminated  or  simply  saved  by  the  computer 
operator  by  depressing  sense  switch  No.  6  on  the  7090/94  console,  or  by 
appropriate  changes  to  the  Data  Definition  (DD)  cards  on  the  S/360.  (See 
subsection  7.4. 3. 2  of  Volume  I).  Using  this  feature,  a  PROGRAM  SAVE  TAPE 
could,  for  example,  be  generated  every  15  minutes  on  long  running  transient 
solutions,  eliminating  the  need  for  repeating  previous  calculations  in  the 
event  of  an  abnormal  run  termination. 

If,  at  the  conclusion  of  a  transient  analysis  run,  the  OUTPUT  SAVE  TAPE 
is  saved,  the  Re-Output  Processor  may  be  used  to  reproduce  lists  and  graphs 
of  any  of  the  circuit  quantities  originally  requested  for  output.  In  addition, 
graphs  of  variables  plotted  against  variables  other  than  time  which  may  not 
have  been  requested  originally  can  be  conveniently  produced. 
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SECTION  IV 


SPARSE  MATRIX  TECHNIQUES 


A  very  common  operation  that  is  performed  in  many  aspects  of  scientific 
digital  computation  is  solution  of  the  system 

A  x  “  b  (118) 

where  A  is  a  coefficient  matrix,  b  is  a  known  vector  of  forcing  functions  and 
x  is  the  vector  of  unknown  quantities  to  be  determined.  Conventional  techniques 
have  of  course  long  been  known.  Gaussian  elimination  or  Gauss-Jordon  methods 
are  standard  in  most  numerical  analysis  texts.  The  need  for  improvement  becomes 
evident  however  when  A  becomes  large.  If  A  is  of  order  n,  these  methods  require 
a  little  more  than  n^  words  of  core  storage.  When  n  =  100,  more  than  10,000 
words  of  core  are  committed  to  this  operation  alone.  Yet  the  larger  A  is,  it 
is  usually  true  that  the  percentage  of  zero  elements,  or  sparseness,  increases. 

A  sparseness  of  90  percent  is  not  unusual  when  n  r  100.  A  sparse  matrix  technique 
then,  can  reasonably  be  defined  as  any  method  that  takes  advantage  of  the 
inherent  sparseness  of  an  often  used  matrix  to  perform  computation  with  some 
significant  saving  of  storage  and  sometimes  even  solution  time. 

The  most  pressing  need  for  some  type  of  sparse  matrix  technique  in  SCEPTRE 
occurs  in  the  DC  algorithm.  The  nature  of  the  mathematical  formulation  is  such 
that  the  A  matrix  corresponding  to  Equation  118  is  equal  in  size  to  the  number 
of  network  resistors  plus  semiconductor  devices.  Experience  has  shown  that 
the  largest  A  matrix  that  could  be  solved  within  200K  bytes  in  the  System/360 
was  about  120  x  120.  This  problem  caused  some  runs  to  fail  d  at  otherwise 
could  have  been  easily  accommodated. 

The  specific  method  that  was  implemented  is  based  upon  triangular  de¬ 
composition  (Ref.  17).  If  A~1  exists,  then  A  can  be  decomposed  as 

A  B  LU  (119) 

where  L  is  a  lower  triangular  matrix  and  U  is  an  upper  triangular  matrix  as 
shown  below: 
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The  matrices  L  and  U  can  be  generated  from  the  original  A  matrix  by  use  of  the 
recursive  relations 


-74- 


0  for  i <  q 


q-i 

2iq  *  aiq  J  Aik  %  with  £iq  “ 

k-1 


u 


Pj 


with  Upj  =  0  for  p>j 


Up j  -  1  for  p  -  j 


Once  L  and  U  have  been  generated,  solution  may  proceed  in  a  straight  forward 
manner.  Since  the  basic  problem  can  be  written  in  the  form  of  Equation  118 
where  vector  X  represents  the  unknowns,  then  substitution  of  Equation  119  in 
118  gives 

L  U  X  =  b  (120) 

Define 

Y  =  U  X  (121) 

which  transforms  Equation  120  to 

L  Y  =  b  (122) 

equation  122  is  in  a  form  suitable  to  solve  for  Y  by  a  simple  front  substitution 

process.  Once  Y  is  known.  Equation  121  can  be  solved  by  back  substitution  to 
arrive  at  the  desired  system  unknown  X.  The  attractiveness  of  this  method  is 
that  it  lends  itself  to  an  indexing  scheme  that  can  work  with  only  nonzero 
elements  of  L  and  non-zero  off  diagonal  elements  of  U.  If  A  is  large  and 
sparse,  a  significant  storage  advantage  can  be  realized. 

The  new  sparse  technique  was  tested  on  two  circuits  that  were  tar  larger  than 
could  have  been  handled  with  the  existing  method  within  200K  bytes  of  core. 

The  first  circuit  was  composed  of  200  resistors  and  semiconductor  junctions 
and  the  resulting  A  matrix  was  well  over  90  percent  sparse.  Accurate  convergence 
was  obtained  in  seven  passes.  The  second  circuit  contained  290  resistors  and 
semiconductor  junctions  and  a  check  of  the  A  matrix  revealed  over  95  percent 
sparseness.  This  circuit  also  converged  accurately  in  seven  passes.  If  this 
same  network  was  solved  with  the  conventional  technique,  approximately  800K 
bytes  of  core  storage  would  have  been  required.  A  series  of  smaller  networks 
were  also  solved  with  the  sparse  technique  and  in  every  case  the  solution 
converged  to  the  same  result  in  the  same  number  of  passes  as  did  the  existing 
method. 
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Since  this  method  is  attractive  for  large  A  only,  it  has  been  programmed 
into  SCEPTRE  such  that  it  becomes  operational  for  networks  containing  70  or 
more  elements  with  DC  analysis  requested.  The  conventional  Gauss-Jordan 
technique  will  continue  to  be  used  for  smaller  networks.  The  sparse  technique 
will  also  be  applied  to  the  transient  portion  of  the  program  for  the  solution 
of  larger  problems.  It  is  clear  that  this  technique  is  a  necessary  prerequisite 
in  eventually  permitting  the  analysis  of  very  large  circuits. 


SECTION  V 


ASSEMBLER  LANGUACE  INPUT/OUTPUT  (I/O) 


SCEPTRE  was  originally  written  in  FORTRAN  IV  on  the  7094  Computing  System 
and  was  later  converted  to  OS/360  FORTRAN  IV  (H-level,  Optimization  =2).  The 
decision  to  program  in  FORTRAN  facilitated  conversion  of  the  program  for  use  on 
other  computing  systems;  but  it  also  resulted  in  increased  CPU  time  requirements 
for  the  ■  tup  phase.  In  an  effort  to  reduce  this  time  requirement  it  was  deter¬ 
mined  .'•.nuL  FORTRAN  I/O  accounted  for  more  than  50  percent  of  the  total  CPU  time 
used,  particularly  for  small  problems.  Since  the  CPU  requirements  for  assembler 
language  I/O  is  at  least  80%  less  than  that  required  by  FORTRAN  I/O  with  FORMAT, 
the  FORTRAN  I/O  was  replaced  with  assembler  language  I/O  with  the  exception  of 
output  directed  to  the  printer.  The  implementation  of  assembler  language  I/O 
into  the  setup  phase  and  the  resulting  improvements  are  discussed  in  the  follow¬ 
ing  paragraphs. 

Three  assembler  language  subroutines  were  written  to  perform  the  functions 
of  the  FORTRAN  7/0  statements  READ,  WRITE,  and  ENDFILE.  These  subroutines 
utilize  GET  and  PUT  assembler  language  instructions  to  perform  the  I/O.  Since 
these  subroutines  transfer  data  to  and  from  core  storage  in  blocks  of  80  bytes 
(characters),  the  logical  record  lengths  ind  the  record  formats  for  the  affected 
data  sets  must  be  80  bytes  and  fixed  block  'FB) ,  respectively.  The  physical 
record  lenj;.  .hs  (i.e.,  BLKSIZES)  can  be  any  multiple  of  80.  The  technique  of 
selecting  the  proper  data  set  for  transferral  of  data  is  via  the  DDNAME  of  the 
data  set.  An  eight  character  DDNAME  for  each  data  set  is  initialized  in  the 
BLOCK  DATA  subroutine.  The  DDNAMEs  used  are  as  shown  in  the  SCEPTRE  Program 
Control  Deck  in  Volume  I,  Figure  81.  If  for  some  reason  the  user  desires 
different  DDNAMEs,  he  need  only  make  the  appropriate  changes  in  BLOCK  DATA 
and  his  SCEPTRE  Program  Control  Deck.  No  reformatting  of  tne  I/O  is  performed 
by  the  assembler  language  subroutines  unless  the  amount  of  data  to  be  trans¬ 
mitted  is  less  than  80  bytes.  In  this  case,  the  remainder  of  the  record  is 
padded  with  blanks. 

The  implementation  of  assembler  language  I/O  did  not  require  an  extensive 
programming  effort  to  reformat  the  I/O  data.  Much  of  the  I/O  of  the  setup 
phase  involves  reading  and  writing  data  where  the  core  storage  and  data  set 
representation  are  both  character  format.  Conversion  to  a  character  format 
and  vice-versa  was  necessary  for  the  data  where  the  core  storage  representation 
was  integer  or  real  constants.  This  includes  circuit  data  stored  in  the  pro¬ 
gram  save  data  set  and  indices  in  the  model  library  data  sets.  The  reformatt¬ 
ing  of  this  data  was  performed  by  FORTRAN  subroutines. 
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An  effort  was  made  to  ensure  that  the  addition  of  the  assembler  language 
I/O  would  not  complicate  the  task  of  converting  the  SCEPTRE/360  System  to 
other  computers.  All  replaced  FORTRAN  I/O  statements  and  their  associated 
FORMAT  statements  were  made  comments  and  left  in  the  programs.  To  assist  in 
the  identification  of  these  statements,  the  character  "X"  was  also  added  in 
column  two  of  the  commented  statements.  Most  of  the  I/O  conversion  effort 
only  requires  that  the  calls  to  the  assembler  language  routines  and  the 
characters  "CX"  from  column  1-2  of  the  FORTRAN  I/O  statements  be  removed. 

Other  statements  to  be  deleted  can  be  easily  identified. 

The  results  of  the  improvement  in  the  execution  speed  of  the  setup  phase 
are  shown  in  table  IV.  The  improvement  was  60  to  75  percent  for  small  problems. 
Although  the  improvement  was  only  30  percent  for  the  large  problem  there  was 
a  significant  reduction  in  total  CPU  time.  As  shown  in  table  VI  the  timing 
tests  were  runs  made  on  an  IBM  System  360  Model  65  computer.  The  finished 
product  will  appear  transparent  to  the  user  in  that  the  I/O  data  itself  will 
retain  its  usual  format. 


TABLE  Iv* 

TIMING  TESTS  OF  THE  SETUP  PHASE 
USING  FORTRAN  I/O  AND  ASSEMBLER  LANGUAGE  I/O 


CPU  Time  Requirements  for 
the  Setup  Phase  (seconds) 

Test  Circuit 

FORTRAN 

I/O 

Assembler 
Language  I/O 

Example  1  from  SCEPTRE 
User's  Manual 

11 

5 

Example  2  from  SCEPTRE 
User's  Manual 

24 

9 

Example  3  from  SCEPTRE 
User's  Manual 

23 

6 

Example  4  from  SCEPTRE 
User's  Manual 

11 

5 

221  Element  Circuit 
(Transient  Solution 
Only) 

90 

62 
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APPENDIX  A 


IMPLICIT  INTEGRATION  IMPLEMENTATION 


A. 1  INTRODUCTION 


Two  basic  factors  are  involved  in  the  amount  of  solution  time  required 
to  complete  any  transient  problem;  the  number  of  steps  or  solution  increments 
needed  and  the  amount  of  computation  required  per  step.  The  first  of  these 
factors  is  controlled  entirely  by  the  type  of  numerical  integration  method 
that  is  used  to  integrate  the  system  of  differential  equations  generated. 

This  system  can  generally  be  expressed  by  the  matrix  equation: 

Y  =  AY  +  BU,  (1) 

0 

where  Y  is  the  vector  of  state  variables,  Y  is  the  vector  of  state  variable 
derivatives  and  both  A  and  B  are  appropriate  coefficient  matrices. 

It  is  well  known  that  the  step  sizes  that  can  be  taken  by  all  explicit 
numerical  integration  methods  are  limited  by  the  largest  system  eigenvalue. 
This  limit  is  often  referred  to  as  the  stability  radius, and  it  differs  from 
method  to  method.  For  example,  if  |\m  |  is  the  absolute  value  of  the  largest 
eigenvalue  of  A  in  Equation  (1),  then  for  most  well  known  explicit  methods 


hm  ^  X  (2) 

T^r 

where  hm  is  the  maximum  step  size  that  may  be  taken.  If  larger  steps  are  taken, 
the  solution  will  begin  to  oscillate.  The  quantity  X  in  Equation  (2)  is  equal 
to  1  for  Euler's  method,  2  for  the  explicit  trapezoidal  method,  2.78  for  Runge 
Kutta  and  6  for  TRAP  2  which  is  currently  an  optional  choice  in  SCEPTRE.  No 
constant  value  for  X  can  be  given  for  XPO,  but  it  too  has  a  stability  limit. 

If  |  \  m  |  is  large  for  any  network,  Equation  (2)  will  force  a  small  step  size 
and  therefore  require  that  many  solution  steps  be  taken  with  an  attendant  rise 
in  computer  solution  time. 

There  has  been  a  recent  trend  in  mathematical  literature  toward  integration 
methods  that  do  not  exhibit  the  limitation  inherent  in  Equation(2).  A  generalized 
form  of  these  methods  can  be  given  as 

P  P 

Vl  ■  I  at  Vi  +h  I  bi  Vi.  (3) 

i=0  i=-l 

where  the  Y  and  h  quantities  remain  as  previously  defined  and  the  a  and  b 
quantities  are  usually,  but  not  necessarily,  constant.  Multistep  methods  are 
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introduced  if  index  Pj<0.  The  simplest  derivative  form  of  Equation  (3)  occurs 
if  coefficients  Sq  ■  1,  b  ^  ■  1  and  all  others  are  set  equal  to  zero.  We  then 
have 


(4) 


which  is  commonly  referred  to  as  the  implicit  or  backward  Euler  technique. 

This  technique  can  be  shown  to  be  unconditionally  stable,  and  has  given  marked 
speed  advantages  over  all  of  the  explicit  methods  in  SCEPTRE  whenever  the 
latter  has  been  affected  by  the  limitation  implied  in  Equation  (2)  A  major 
disadvantage  of  this  particular  method  lies  in  its  inaccuracy  for  some  types 
of  problems.  Implicit  or  not,  it  is  still  a  first  order  method  and  tests 
have  indicated  that  both  XPO  and  RUK  can  be  considerably  more  accurate.  The 
next  well  known  variant  of  Equation(3)  is  set  up  if  coefficients  a^  =  1, 
b  ^  =  0.5  and  all  others  are  set  to  zero.  This  choice  sets  up 


v  =  y  +  —  y  +  —  Y 

n+l  n  2  n  2  n+1,  (5) 

a  method  called  implicit  trapezoidal  integration.  This  method  is  considerably 
more  accurate  than  the  first  order  variety,  but  has  lost  the  quality  of 
unconditional  stability.  Under  some  computation  circumstances  oscillations 
can  occur  and  render  the  analysis  suspect.  Therefore,  after  considerable 
testing  and  consultation  with  the  Air  Force  Weapons  Laboratory,  it  was  decided 
to  implement  a  multistep  (Ref.  8,9)  version  of  Equation  (3)  that  automatically 
chooses  Is  PS  6.  A  good  deal  of  testing  has  been  performed  on  this  implementation 
and  the  resultant  conclusions  are  given  in  the  following  paragraphs. 

A. 2  MANUAL  EXAMPLE  1 

The  first  practical  network  to  be  run  on  SCEPTRE  with  implicit  integration 
was  example  1  in  the  SCEPTRE  manual  Vol.  I,  Section  4.1.  This  particular  net¬ 
work  contains  a  significant  spread  in  eigenvalues  due  to  the  small  transistor 
capacitors  and  the  larger  load  capacitor.  A  master  run  with  two  associated 
reruns  was  made  and  a  summary  of  the  number  of  integration  steps  and  passes 
required  by  both  the  original  explicit  method  and  the  new  implicit  method  is 
given  in  Table  IV.  Since  this  network  did  exhibit  a  significant  spread  in 
eigenvalues,  the  implicit  results  were  much  faster.  No  practical  difference 
in  accuracy  was  noted  in  any  part  of  the  runs. 


Note:  Large  Eigenvalues  can  inhibit  the  step  size  and  small  Eigenvalues  can 
make  a  large  problem  duration  necessary.  Therefore  it  is  the  spread 
in  Eigenvalues  that  often  causes  difficulty  for  explicit  integration 
methods . 
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TABLE  A-l 


COMPARISON  OF  EXPLICIT  AND  IMPLICIT  INTEGRATION 
ON  EXAMPLE  1 


Explicit 

Implicit 

Passes 

Steps 

Passes 

Steps 

Master 

6441 

3049 

291 

126 

First  Rerun 

6065 

2867 

480 

155 

Second  Rerun 

5972 

2832 

255 

112 

A. 3  A  PRACTICAL  ILLUSTRATION 


The  next  network  to  be  discussed  delves  a  bit  more  deeply  into  the 
practical  situation  that  gives  implicit  integration  an  advantage  over 
explicit  integration.  Consider  the  circuit  in  Figure  19  in  which  a  two- 
stage  transistor  circuit  drives  a  large  capacitive  load.  Let  the  equivalent 
circuit  component  data  be  such  that  the  emitter  and  collector  capacitances 
are 


CE  = 

Cte 

+  e  t  j 

nee 

(6) 

cc  = 

Ctc 

+  0,  T  J 
i  s  c 

(7) 

where 

C  * 
te 

= 

emitter  transition  capacitance  = 

3  Pf 

Ctc* 

= 

collector  transition  capacitance 

=  5  pf 

0n 

= 

30  V-1 

0i 

= 

35  V_1 

T 

e 

= 

0.2  ns  -*  f Q  =  800  me 

T 

s 

= 

5  ns 

J 

e 

= 

forward  emitter  junction  current 

J 

c 

= 

forward  collector  junction  current 

if 

The  transition  capacitances  are  assumed  to  be 
does  not  affect  the  point  of  the  discussion. 

constant  here. 

This  simplification 
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+  10  V 


Figure  19.  A  Practical  Illustration 

If  the  input  signal  contains  high  frequency  components,  or  if  the  performance 
of  the  individual  transistor  stages  are  of  interest,  the  user  cannot  "cheat" 
on  the  equivalent  capacitor  values.  The  available  component  data  must  be  used 
in  Equations  (6)  and  (7)  even  if  these  lead  to  very  small  values  of  capacitance. 
On  the  other  hand  the  values  used  for  Rl3  and  lead  to  an  output  time  constant 
of  2  ms  which  in  turn  requires  a  problem  duration  of  about  10  ms  if  the  entire 
transient  is  to  be  observed.  The  familiar  eigenvalue  spread  has  been  created. 
When  this  problem  was  solved  with  explicit  integration  (SPO) ,  10,000  passes 
were  required  to  get  to  a  problem  time  of  58,000  ns  -  far  short  of  the  required 
duration.  Solution  with  the  new  implicit  method  covered  the  entire  transient 
in  520  passes.  With  this  data  one  can  claim  an  improvement  factor  of  3260. 

This  is  based  on  a  factor  of  19.2  fewer  passes  to  achieve  a  problem  duration 
of  another  factor  of  170  greater.  Entire  papers  have  been  written  that  were 
based  on  statistics  like  this. 
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The  same  network  can  be  used  to  illustrate  that  there  is  still  a  future 
for  explicit  integration.  If  components  Rl3  and  are  removed  from  the 
network,  so  is  the  large  time  constant  and  any  need  to  run  the  problem  for 
anything  like  10  ms.  In  fact,  100  ns  is  more  than  sufficient  to  cover  the 
entire  transient  period.  What  has  been  done  in  effect  is  to  remove  the 
large  spread  in  eigenvalues.  This  change  in  computational  circumstances 
also  changes  the  relative  efficiencies  of  the  integration  methods.  Explicit 
integration  now  completes  the  problem  in  297  passes  while  implicit  requires 
409  passes.  Without  an  eigenvalue  spread,  the  implicit  method  has  lost  its 
large  speed  advantage  and  has  become  about  40  percent  slower. 

The  complete  input  listing  of  the  circuit  of  Figure  19  including  the 
model  description  is  given  on  the  next  page.  It  is  interesting  to  note  that 
the  Rerun  mode  may  be  used  to  easily  perform  comparisons  between  implicit 
integration  and  explicit  XP0. 

A. 4  EIGENVALUE  SPREADS 


A  series  of  runs  and  reruns  were  made  in  order  to  gather  comparative  data 
on  the  performance  of  XPO  and  the  implicit  method  as  a  function  of  selected 
eigenvalue  spreads.  This  information  is  summarized  in  Table  A-II.  With  a  small 
spread  of  ten,  it  is  seen  that  implicit  takes  more  solution  passes  than  does 
explicit  -  hence  an  improvement  factor  of  less  than  unity.  As  the  spreads 
of  eigenvalues  are  progressively  increased,  implicit  integration  becomes 
comparatively  more  attractive, and  the  trend  would  become  more  marked  if  the 
spread  sequence  had  been  continued  on. 

The  forcing  function  used  in  all  runs  was  a  1  V  step  function  in  order 
to  avoid  complicating  effects  from  that  direction.  If  more  complex  forcing 
functions  were  used,  it  is  likely  that  the  improvement  factors  that  are  listed 
in  Table  III  would  be  modified  and  in  some  cases  quite  significantly.  Since 
this  is  true,  the  listed  factors  should  be  considered  to  be  valid  only  as  an 
approximate  guide  that  illustrates  the  potential  of  implicit  integration 
under  favorable  cirrumstances . 


TABLE  A-II 


EXPLICIT-IMPLICIT  COMPARISON  AT 
SELECTED  EIGENVALUE  SPREADS 


Eigenvalue 

Spread 

(approximate) 

Explicit 

(XPO) 

Passes 

Implicit 

Passes 

Improvement 

Factor 

10 

100 

1000 

10000 

100000 

10  o 

10  7 

108 

113 

237 

703 

1807 

6664 

16578 

*20000 

*20000 

146 

171 

180 

190 

184 

185 

193 

189 

0.77 

1.4 

3.9 

9.5 

36.0 

90.0 

238.0 

996.0 

*The  last  two  runs  were  not  carried  through  the  complete  transient 
solution,  but  were  terminated  by  the  program  pass  limit.  The 
improvement  factors  for  these  runs  include  an  extrapolation  term 
that  is  based  on  the  portion  of  the  run  that  was  not  completed. 
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MODEL  DESCRIPTION 
MODEL  ZZ  (TEMP)  (B-C-E) 

ELEMENTS 

RB, B-X=.3 

RC, C-Y=.015 
CE,X-E=X1(3.+6.*JE) 
CC,X-Y=X2(5.+175.*JC) 

JE ,X-E=DIODEQ(l .D-7,30.) 

JC  ,X-Y=DIODEQ(l .D-7,35. ) 

JX, Y-X=.98*JE 

JY, E-X=.1*JC 
OUTPUTS 
VCE,VCC,PLOT 
CIRCUIT  DESCRIPTION 

TWO  STAGE  CIRCUIT  WITH  LOAD 

ELEMENTS 

El , 1 -2=T ABLE1 

RBI  ,2-3=2 

T1 ,3-4-5=M0DEL  ZZ 

E2, 1-4=10 

RL2,4-7=.5 

RE ,  1 -5= . 2 

RB2 ,5-6= . 5 

T2, 6-7-1 =M0DEL  ZZ 

RL3, 7-9=1 

CL,9-1=2E6 

OUTPUTS 

VCL  ,XSTPSZ  ,PL0T 

RUN  CONTROLS 

RUN  INITIAL  CONDITIONS 

INTEGRATION  ROUTINE= IMPL ICIT 

STOP  TIME=1E7 

FUNCTIONS 

TABLE1 

0,0,  50,5,  100,5 

RERUN  DESCRIPTION 

RUN  CONTROLS 

INTEGRATION  ROUTINE=XPO 

MAXIMUM  INTEGRATION  PASSES=10000 

END 
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A.  5  JACOBIAN  CONSTRUCTION 


The  need  for  a  Jacobian,  or  matrix  of  partial  derivatives,  arises  in  the 
derivation  of  the  corrector  iteration  of  the  implicit  integration  method  that 
has  been  implemented.  This  need  is  common  to  all  implicit  methods,  both 
single  step  and  multistep.  An  m  dimensional  generalized  form  of  the  Jacobian 
appears  in  Figure  20.  Two  different  philosophies  exist  concerning  the  most 
efficient  method  of  forming  this  matrix.  The  first  is  based  on  a  symbolic 
approach  in  which  each  element  of  the  Jacobian  is  set  up  in  terms  of  the  literal 
R,  C,  L  and  G  components  that  constitute  a  given  problem.  Its  advantage  lies 
in  the  fact  that  the  matrix  can  be  constructed  just  once  and  updated  as  often 
as  required  depending  on  the  progress  of  the  solution.  Furthermore,  the  up¬ 
dating  process  is  just  a  simple  matter  of  data  insertion  without  any  need  to 
recompute  the  system  derivatives.  When  large  problems  are  encountered,  sparse 
matrix  techniques  can  be  applied  without  difficulty.  The  disadvantage  is 
that  some  combinations  of  topologies,  dependencies,  and  element  values  are 
encountered  such  that  the  Jacobian  that  is  constructed  this  way  is  significantly 
in  error.  In  these  cases  the  average  solution  step  size  that  can  be  taken  is 
sometimes  significantly  reduced  and  inefficient  operation  can  result. 
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Figure  20.  The  General  Jacobian 


The  second  method  is  purely  numerical  in  nature  and  can  be  referred  to 
as  "numerical  differencing".  It  produces  an  approximation  to  the  desired 
partial  derivatives  as  indicated  in  Equation  (8).  These  approximations  are 
made  by  making  m 


0 


Y  ) 
m,t 


(8) 


additional  derivative  computations  each  time  the  Jacobian  is  to  be  updated. 
The  advantage  of  this  approach  is  that  it  is  universal  in  that  theoretically 
it  should  produce  a  good  approximation  to  the  true  Jacobian  for  any  transient 
problem.  Its  principal  disadvantage  is  that  it  requires  additional  solution 
passes.  The  number  of  additional  solution  passes  will  be  equal  to  mX  where  X 
is  the  number  of  Jacobian  evaluations  required  during  the  course  of  a  given 
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problem.  A  second  disadvantage  is  Chat  this  method  does  not  readily  lend 
Itself  to  the  use  of  the  L-U  decomposition  sparse  matrix  technique  that  is 
used  in  the  program. 

The  program  will  automatically  choose  between  the  two  methods  of  Jacobian 
construction.  The  procedure  that  is  followed  is  to  categorize  all  the  transient 
problems  that  are  to  be  solved  with  implicit  integration  into  three  groups 
according  to  the  number  of  state  variables  contained,  m. 

Group  1  -  m  ^ 10 

Group  2  -  10 <  m  ^50 

Group  3  -  m  >50 

The  program  default  is  set  to  use  the  numerical  approach  for  Group  1 
problems.  The  reason  for  this  is  th  *;  both  of  the  disadvantages  of  this 
approach  (m  X  extra  derivative  evaluations  and  additional  core  requirements) 
are  greatly  minimized  for  small  problems.  Exactly  the  convt  rse  is  true  for 
Group  3  problems  so  the  default  choice  then  becomes  the  symbolic  approach. 

When  the  medium  size  problems  represented  by  the  Group  2  classification  are 
encountered,  a  series  of  checks  will  be  made.  These  checks  will  examine  the 
problem  for  topological  conditions  that  could  lead  to  error  as  the 

presence  of  user  supplied  differential  equations.  If  any  of  >  *s  are 

positive  the  numerical  approach  will  be  taken.  If  not,  symbo  ‘  .  action 

of  the  Jacobian  will  be  used. 

A  method  has  also  been  provided  to  allow  the  user  to  specify  '.,.e  type  of 
construction  to  be  used  which  will  bypass  the  default  mechanism.  If  either 
the  numerical  or  symbolic  approach  is  desired  the  respective  entry  under  RUN 
CONTROLS  is 

USE  DIFFERENCED  JACOBIAN 

USE  SYMBOLIC  JACOBIAN 
A.  6  PRESET  STEP  SIZE  CONTROLS 

The  explicit  integration  routines  in  SCEPTRE  have  preset  quantities 
that  control  the  maximum,  minimum,  and  starting  step  sizes.  These  are 
identical  for  all  of  the  explicit  methods  and  are  as  follows: 

MINIMUM  STEP  SIZE  =  1  x  10~5  (STOP  TIME) 

MAXIMUM  STEP  SIZE  =  2  x  10_2  (STOP  TIME) 

STARTING  STEP  SIZE  =  1  x  10"3  (STOP  TIME) 

The  preset  quantities  have  been  found  to  be  inappropriate  for  implicit 
integration  and  have  been  replaced  with  the  following: 

MINIMUM  STEP  SIZE  =  1  x  10"14  (STOP  TIME) 

MAXIMUM  STEP  SIZE  -  2  x  10_2  (STOP  TIME) 

STARTING  STEP  SIZE  -lx  10"8  (STOP  TIME) 
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Any  of  these  may  be  replaced  with  a  different  constant  at  the  discretion  of 
the  user  with  a  simple  entry  under  RUN  CONTROLS.  This  format  is 

MINIMUM  STEP  SIZE  -  number 
MAXIMUM  STEP  SIZE  -  number 
STARTING  STEP  SIZE  “  number 

The  error  control  differs  most  markedly  from  that  used  for  the  explicit 
methods  in  that  the  estimated  error  is  not  computed  and  checked  for  each 
individual  differential  equation;  but  is  computed  and  checked  for  all  of 
the  equations  together  just  once  each  step.  If  a  given  problem  consists  of 
m  state  variables,  the  following  relation  is  computed  at  the  end  of  each 
integration  step: 


2 

E. 

E 

1 

+  . . .  + 

m 

lYi' 

|Y  | 

1  m' 

_ 

- 

where  the  Ei  terms  designate  local  error  estimates  for  each  differential 
equation  and  the  Yi  are  the  state  variables  with  0<i£m.  The  step  is 
accepted  if 

DSg  *MINIMUM  ABSOLUTE  ERROR 

where  g  is  a  constant  that  depends  on  the  current  order  of  the  integration 
method  and  MINIMUM  ABSOLUTE  ERROR  has  a  preset  value  of  0.001.  The  nomenclature 
notwithstanding,  this  is  a  relative  error  control  since,  the  magnitude  of  the 
state  variables  enter  the  computation. 

Some  testing  has  been  performed  to  determine  the  effect  of  other  values 
for  MINIMUM  ABSOLUTE  ERROR.  No  significant  improvement  was  found  in  the 
inevitable  speed-accuracy  tradeoff.  What  was  found  however,  was  that  larger 
values  which  lead  to  significant  speed  improvements  did  so  only  at  the  cost 
of  unacceptable  degradation  in  accuracy.  The  user  is  free  to  experiment 
along  these  lines  by  entering  under  RUN  CONTROLS,  MINIMUM  ABSOLUTE  ERROR  ■ 
number.  The  remaining  error  criteria  (MAXIMUM  ABSOLUTE  ERROR,  MAXIMUM 
RELATIVE  ERROR  and  MINIMUM  RELATIVE  ERROR)  will  have  no  effect  on  the  implicit 
integration  method. 
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APPENDIX  3 
B  MATRIX  DERIVATION 


The  derivation  of  the  general  B  matrix  that  expresses  link  voltages  in 
terms  of  tree  branch  voltages  and  tree  branch  currents  in  terms  of  link 
currents  is  as  follows: 


Two  fundamental  incidence  matrices  that  arise  from  network  topology 
theory  will  be  called  the  Q  and  T  matrices  here.  The  fundamental  cut  set 
matrix,  Q  =  t^ijl  *s  a  matrix  containing  (n-1)  rows  and  b  columns  for  a  net¬ 
work  containing  n  nodes  and  b  elements,  where: 


qij 

qij 

qU 


+1  if  the  i1"*1  fundamental  cut  set  direction  coincides  with  the 
reference  direction  of  the  jth  element 

-1  if  the  i^  fundamental  cut  set  direction  is  in  opposition  to 
the  reference  direction  of  the  j element 

0  if  the  i1"*1  fundamental  cut  set  does  not  include  the  j1"*1  element 


If  the  elements  are  properly  ordered,  it  is  always  true  that 


Q  -  [-BT  u] 

where  the  columns  of  the  unit  matrix  U  correspond  to  the  tree  branch  elements. 


The  fundamental  circuit  matrix,  T  =  [tij]  is  a  matrix  containing  m  rows 
and  b  columns  for  a  network  containing  m  independent  loops  and  b  elements, 
where. 


+1  if  the  ith  independent  loop  direction  coincides  with  the 
reference  direction  of  the  element 


tij  =  -1  if  the  ita  independent  loop  direction  opposes  the  reference 
direction  of  the  element 

t^j  =  0  if  the  i^  independent  loop  does  not  include  the  jth  element 


If  the  elements  are  properly  ordered,  it  is  always  true  that 
T  =  [U  B] 


where  the  columns  of  the  unit  matrix  U  correspond  to  the  network  links. 
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Since  it  is  always  true  that 


QIb  «  0,  T  Vb  -  0 


direct  substitution  yields 


Expansion  of  these  relations  gets 
ITB  -  bt  il  and  VL  >  -B  VTB 

so  that  the  tree  branch  currents  may  be  expressed  in  terms  of  the  link  currents 
and  the  link  voltages  may  be  expressed  in  terms  of  the  tree  branch  voltages 
through  the  B  matrix. 


APPENDIX  C 

DIODE  REPRESENTATION  IN  THE  INITIAL-CONDITIONS  PROGRAM 


The  current  In  a  diode  or  transistor  junction  at  a  point  on  its  voltage 
vs.  current  (V/l)  characteristic  may  be  separated  into  two  components  as 
J  ■  GjVyK}.  Consider  the  typical  diode  curve  below. 


Angle  Oj  =  angle  a*  is  enclosed  by  the  slope  of  the  diode  characteristic 
at  any  point  and  the  horizontal  at  that  point.  E'  is  an  offset  voltage  that 
marks  the  intersection  of  a  continuation  of  the  slope  Gj  and  the  line  J  •  0. 

From  the  figure: 


tan  a. 


tan  a„  *  G. 


Vj  -  E’ 


or  J 


GJ  VJ  -  GJ  E’ 


or  J  =  Gj  Vj  +  Q  if  Q  is  defined  as  -G^E' 
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APPENDIX  D 


BASIC  NEWTON-RAPHSON  METHOD 

Given  a  single  algebraic  or  transcendental  equation  of  the  form  F(x)  ■  0, 
that  is  single  valued  and  differentiable  in  the  domain  of  interest,  a  Newton- 
Raphson  procedure  can  be  constructed  using 

F(x)  +  Ax  =  0 

a  x 

An  initial  x  =■  xq  may  be  assumed  and  F(xq), 


9F(xq) 


determined. 


Then  A  x 


SF(x0)  “d  X1  ‘  xo  +  A* 
9x 


The  procedure  is  repetitive  until 
F  (xn+1)  -F(xnS)  !  <  z 


where  z  is  some  specified  convergence  criteria.  When  this  last  relation  is 
satisfied,  the  process  is  said  to  have  converged.  Extension  to  systems  of 
equations  adds  no  complications. 
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APPENDIX  E 

EQUIVALENCE  OF  EQUATIONS  (72)  AND  (72') 


To  show  the  equivalence  of  equations  (72)  and  (72'),  it  will  be 
sufficient  to  show  that 


"V9 

F1  <V9,  V2.  v5> 

-1 

-b97  e7 

V2 

-[z]"1 

F2  (V9,  V2,  V 

-tz] 

-B27  E 7 

r  m  _  i 

V5 

f3  (V9j  v2,  vs) 

[b95  T  +  B05  TqJ  q9  +  B85  TJ8 

or 


V9 

-b97  e7 

"F1  (V9,V2,V5  >‘ 

V2 

*B27  E7 

+ 

F2  (V9,V2,V5  > 

Vr 

'  T  _  T1  „  T, 

_  5_ 

B95  +B05  °J  ^9  +  B85  J8_ 

F3  <V9,V2,  V5> 

* 

or 


V9 

-B97  E7 

’F1 

(V9.V2, 

V" 

[Z] 

V2 

=3 

-B27  E7 

r 

+ 

F2 

<V9,V2, 

V 

_V5 

[B95T  +  B05Tq]  «9  +  B85Tj8 

F3 

<V9,V2, 

V 

The  left  side  of  the  equation  expands  into 


which  shows  the  equivalence. 


-93- 


APPENDIX  F 


ADJOINT  NETWORK  FOR  SENSITIVITY 


A  number  of  papers  (Ref.  10,  11)  have  been  written  on  determing  the  un- 
normalized  sensitivity  (i.e.,  partial  derivatives)  of  network  functions  with 
respect  to  network  elements  by  an  adjoint  network.  To  find  the  sensitivity 
of  any  network  function  by  this  method,  two  analyses  are  required,  one  of  the 
given  network  and  the  other  of  the  related  adjoint  network.  This  appendix 
describes  the  generation  of  the  adjoint  network  for  DC  calculations  in  SCEPTRE. 
A  detailed  derivation  of  the  adjoint  network  for  sensitivity  is  given  in 
References  (10)  and  (11).  Due  to  the  sign  conventions  followed  in  SCEPTRE, 
the  terms  given  in  tables  F-I  and  F-II  differ,  in  sign,  from  those  given  by 
Director  and  Rohrer.  (Ref.  10,  11)  Table  F-I  describes  the  relation  between  a 
given  network  and  its  related  adjoint  network.  It  also  describes  the  allowed 
dependent  variables  (network  functions)  and  independent  variables  (network 
elements).  Table  F-II  describes  the  forcing  functions  for  the  adjoint  network 
and  the  sensitivities.  The  superscript  "a"  in  the  tables  indicate  the  currents 
and  voltages  in  the  adjoint  network. 
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TABLE  F-I 

RELATION  BETWEEN  NETWORK  AND  ADJOINT  NETWORK 
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Superscript  'a"  indicates  value  frcn  Adjoint  network  solution 
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Secondary  current  sources  in  the  given  network  introduce  voltage 
dependent  voltage  sources  in  the  adjoint  network  (shown  in  Table  f-I) . 

This  requires  modification  of  the  Jacobian.  Derivation  of  the  modified 
Jacobian  follows.  Let  us  classify  the  voltage  sources,  arising  in  the 
adjoint  calculation  due  to  secondary  sources  (J0's)  in  the  main  circuit, 
as  Ejj's. 

The  B  matrix  from  the  L  tree  which  includes  the  ED's  is  as  shown  below. 


Class  C^ 

Class 

Class  Lg 

Class  E^ 

Class  Ed 

Class  C^ 

B14 

B15 

B16 

B17 

bid 

Class  R2 

0 

B25 

B26 

B27 

B2D 

Class 

0 

0 

B36 

B37 

B3D 

Class  Jg 

CO 

B85 

B86 

B87 

B8D 

Class  Jg 

B94 

B95 

B96 

B97 

B9D 

Class  J^ 

B04 

B05 

B06 

B07 

b0d 

The  following  equations  arise  from  the  above  B  matrix  if  vectors  Vg, 
V^,  1^,  1^  and  submatrix  B^  are  assumed  to  be  zero.  These  assumptions 

are  based  on  the  known  final  values  of  V, ,  V_,  I.  and  I,  for  the  initial 

6  3  4  1 

condition  problem  and  on  the  absence  of  any  current-source  and  capacitor 
cut  sets. 

*5  -  -  B85J8 ' 4J9  -  B!)5J0  '  °  (1) 


V2  +  B25V5  +  B27E7  +  B2DED 


V9  +  B95V5  +  B97E7  +  B9DED 


To  express  Equation  (1),  (2)  and  (3)  in  terms  of  V,.,  V- ,  VQ,  Jfl 

?  1.1 _ r.  1 1 _ _ • _ _ _ a.  j _ _  _  _ _» _ i  j  ^  y  O 


and  E^  the  following  equations  are  required. 


*5  -  G55V5 


*2  *  G22V2 
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J9  "  G99V9  +  <>9 


(6) 

(7) 


For  the  definition  of  G^_,  G^t  q,  Ggg  and  Qg,  see  subsection  2.4.1. 
Two  more  equations  required  for  complete  derivations  are  (8)  and  (9)  below. 

E  -  a  *VJ  (8) 

D  0 

where  a*  is  the  diagonal  matrix  whose  (i,i)  element  is  the  non-zero  entry 
fch 

of  the  i  row  of  (Equation  7),  i.e.  it  is  the  coefficient  of  the  Jg  on 
which  the  i*"^  Jq  depends. 

The  absence  of  any  current-source  cut  set  in  the  given  circuit  makes 
Bqd  equal  to  zero.  This  gives... 


VJ0  +  B05V5  +  B07E7 


(9) 


From  Equations  (1)  through  (9),  we  get... 


°55V5  '  »25G22V2  '  B85J8  '  'B95  +  BS5  ° 1  'G99V9  +  V  ’  0 
V2  +  'B25  -  B2D  °*305JV5  +  'B27  "  B2D  °‘B07)E7  '  0 
V9  +  [B95  ’  B9D  “*B05,V5  +  [B97  '  B9D  °‘B07)E7  '  0 

Following  through  a  derivation  similar  to  that  done  in  subsection  2.4.1 
we  get  Equation  (10).  Equation  (11)  is  the  original  SCEPTRE  derivation. 


V5(n+1) 

°55  -B25G22  'B95+BJ>50'099 

1 

1B95+BS5QlVBLJ8) 

V2 (n+1) 

- 

tB25'B2D  °  ‘B05>  1  ° 

>-B27+B2D  “*B07>E7 

V9(n+1) 

tB95"B9D  °  B05]  0  1 

[-B97+B9D°‘B071E7 

(10) 
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Comparison  of  Equations  (10)  and  (11),  shows  that  terms  get  added  to  the 
(2,1)  and  (3,1)  terms  of  the  Jacobian  and  to  the  forcing  function  column 
vector  of  equation  (10),  while  the  terms  originally  involving  Ci  are  deleted. 

Therefore,  when  an  adjoint  run  is  requested,  and  secondary  sources  are 
present  in  the  given  circuit,  the  iteration  given  by  Equation  '10)  is  used 
instead  of  the  one  given  by  (11)  . 
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